2.2 另一个现实例子
数据包中的co2
>
m1.co2=arima(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12)) > m1.co2
Call:
arima(x = co2, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
Coefficients:
ma1 sma1 -0.5792 -0.8206 s.e. 0.0791 0.1137
sigma^2 estimated as 0.5446: log likelihood = -139.54, aic = 283.08
> detectAO(m1.co2) [1] \> detectIO(m1.co2) [,1] ind 57.000000 lambda1 3.752715
拟合含有新息异常的模型
>
m4.co2=arimax(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),io=c(57)) > m4.co2
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-57 -0.5925 -0.8274 2.6770 s.e. 0.0775 0.1016 0.7246
sigma^2 estimated as 0.4869: log likelihood = -133.08, aic = 272.16
模型显示AIC相比之前模型一更小了。而且IO效应的P 值=2.677/0.7246是显著的.
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中可以看出两样本在滞后2期存在明显的相关性。
3.2 奶产量与对数化发电量的伪相关
data(milk)
data(electricity)
milk.electricity=ts.intersect(milk,log(electricity))# intersect函数将多个时间序列合并在一个容器中。
ccf(as.numeric(milk.electricity[,1]),as.numeric(milk.electricity[,2]),
main='milk & electricity',ylab='CCF')
两者相关性似乎非常的强，但实际上这是因为他们的各自存在很强的自相关性。
4 预白化与随机回归
对于具有强自相关的数据而言，很难评估两个过程之前是否存在依赖关系，因而，宜将x和y之间的线性关系关联从其各自相关关系中剥离出来。预白化正是为了达到此目的的一个有效工具。
4.1 牛奶与电量的CCF预白化校正
> data(milk) >
me.dif=ts.intersect(diff(diff(milk,12)),diff(diff(log(electricity),12))) >
prewhiten(as.vector(me.dif[,1]),as.vector(me.dif[,2]),ylab='CCf')
再次分析两者的相关性，此时除了时滞-3具有边缘显著外，其他地方没有一个相关系数是显著的。幌动防震 这给出的35个样本互相关系娄中大约会出现 1.75=35x0.05个虚假警报，即这个-3系数的显著可能就是一个虚假的信息。因此，牛奶与耗电量序列实际上是基本不相关的。从而认为之前在原始数据序列中发现的强互相关是伪相关的。
4.2 Log（销售量）与价格数据的相关性分析 4.2.1 预白化处理
plot(bluebird,yax.flip=T)#画两者的时间序列对比图
预白化处理
prewhiten(y=diff(bluebird)[,1],x=diff(bluebird)[,2],ylab='ccf')
从CCF图可以看出两者之间只在时滞0处是显著的。即价格与销售量之间存在着很强的同期负相关关系。即当期提高价格将导致销售量的当期下降。
4.2.2 一般线性回归分析
> sales=bluebird[,1] > price=bluebird[,2]
> chip.m1=lm(sales~price) > summary(chip.m1)
Call:
lm(formula = sales ~ price)
Residuals:
Min 1Q Median 3Q Max