时间序列回归模型R实现.docx
- 文档编号:30194340
- 上传时间:2023-08-07
- 格式:DOCX
- 页数:49
- 大小:319.43KB
时间序列回归模型R实现.docx
《时间序列回归模型R实现.docx》由会员分享,可在线阅读,更多相关《时间序列回归模型R实现.docx(49页珍藏版)》请在冰豆网上搜索。
时间序列回归模型R实现
**
时间序列回归模型
干预分析
概念及模型
Box和Tiao引入的干预分析提供了对于干预影响时间序列的效果进行评估的一个框架,假设干预是可以通过时间序列的均值函数或者趋势而对过程施加影响,干预可以自然产生也可以人为施加的,如国家的宏观调控等。
其模型可以如下表示:
其中mt代表均值的变化,Nt是ARIMA过程。
干预的分类
阶梯响应干预
脉冲响应干预
**
干预的实例分析
模型初探
对数化航空客运里程的干预模型的估计
>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〕
系数。
**
拟合带有干预信息的模型
函数:
arimax(x,order=c(0,0,0),seasonal=list(order=c(0,0,0),period=NA),xreg=NULL,=TRUE,=TRUE,fixed=
NULL,init=NULL,method=c("CSS-ML","ML","CSS"),n.cond,=
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:
ma1
sma1
Dec96
Jan97
Dec02
I911-MA0
s.e.
sigma^2estimatedas0.0006721:
loglikelihood=219.99,aic=画图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/11Effects')abline(h=0)
从上图可以看出在2003年底后,911事件的影响效应才平息,航班客运量恢复了正常。
**
异常值在时间序列中异常有两种,可加异常和新息异常,分别记AO和IO。
异常值例如
模拟数据
模拟一般的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)y
[9]
[17]
[25]
**
[57]
[65]
[73]
[89]
[97]
>y[10]<-10
模型初步判断
>acf(y)
**
>pacf(y)
eacf(y)
AR/MA
0xxooooooooo
ooo
1ooooooooooo
ooo
2ooooooooooo
ooo
3oxooooooooo
o
o
o
4oxooooooooo
o
o
o
**
5xxooooooooo
ooo
6xoooooooooo
o
oo
7oxooooooooo
o
oo
从三个的结果来看,可以初步分析
y是AR〔1〕模型
对模型时行拟合
m1=arima(y,order=c(1,0,0))
m1
Call:
arima(x=y,order=c(1,0,0))
Coefficients:
ar1
intercept
s.e.
对模拟模型进行异常值探测
>detectAO(m1)[,1][,2][,3]
ind
lambda2
**
>detectAO(m1,robust=F)[,1]
ind
lambda2>detectIO(m1)
[,1]
[,2]
ind
lambda1
AO
探测结果认为第
9,10,11.可能出现异常值。
IO
探测认为第
10,11
可能出现了异常
值。
由于检验统计量的最大取值出现在
10且AO〉IO,所以更认为出现异常值在第
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:
ar1interceptAO
**
s.e.
sigma^2estimatedas1.059:
loglikelihood=-145.29,
aic=
>detectAO(m2)[1]"NoAOdetected"
>detectIO(m2)[1]"NoIOdetected"
比拟有无异常值的两模型
再次进行异常值探测时,没有发现异常值,验证最初序列异常出现在10的猜想
比照模型1和2的拟合效果
>tsdiag(m2)>tsdiag(m1)
**
虽然模型二的残差通过引入异常值后正太性是显性的,但是其acf和P值结果显示引入MA〔1〕是必要的。
重新拟适宜当模型
m3=arima(y,order=c(1,0,1),xreg=data.frame(AO=seq(y)==10))
detectAO(m3)"NoAOdetected">detectIO(m3)"NoIOdetected">tsdiag(m3)>m3Call:
arima(x=y,order=c(1,0,1),xreg=data.frame(AO=seq(y)==10))
**
Coefficients:
ar1
ma1
intercept
AO
s.e.
sigma^2estimatedas0.793:
loglikelihood=-131.16,aic=模型的拟合效果是显著提高。
Acf和P值检验也一步通过。
>plot(y,type='b')arrows(40,7,11,9.8,length=0.8,angle=30)
**另一个现实例子
数据包中的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:
ma1sma1
s.e.
**sigma^2estimatedas0.5446:
loglikelihood=-139.54,aic=
detectAO(m1.co2)[1]"NoAOdetected"detectIO(m1.co2)
[,1]ind
lambda1拟合含有新息异常的模型
>
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:
ma1sma1IO-57
**
s.e.
sigma^2estimatedas0.4869:
loglikelihood=-133.08,
aic=
模型显示
AIC
相比之前模型一更小了。
而且
IO
效应的
P值
是显著的
.
伪相关在时间序列中引入协变量,如非洲牧草产量通常与某些气候指标密切相关,在这种发问下在通过在时间序列模型中纳入相关的协变量,将有助于更好的了解根底过程以及得到更为准确的预测。
模拟数据
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期存在明显的相关性。
奶产量与对数化发电量的伪相关
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')两者相关性似乎非常的强,但实际上这是因为他们的各自存在很强的自相关性。
预白化与随机回归对于具有强自相关的数据而言,很难评估两个过程之前是否存在依赖关系,因而,宜将x和y之间的线性关系关联从其各自相关关系中剥离出来。
预白化正是为了到达此目的的一个有效工具。
牛奶与电量的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个样本互相关系娄中大约会出现个虚假警报,即这个-3系数的显著可能就是一个虚假的信息。
因此,牛奶与耗电量序列实际上是基本不相关的。
从而认为之前在原始数据序列中发现的强互相关是伪相关的。
Log〔销售量〕与价格数据的相关性分析
预白化处理
plot(bluebird,yax.flip=T)#画两者的时间序列比照图
预白化处理
**
prewhiten(y=diff(bluebird)[,1],x=diff(bluebird)[,2],ylab='ccf')
从CCF图可以看出两者之间只在时滞0处是显著的。
即价格与销售量之间存在着很强的同期负相关关系。
即当期提高价格将导致销售量的当期下降。
一般线性回归分析
sales=bluebird[,1]price=bluebird[,2]chip.m1=lm(sales~price)summary(chip.m1)
Call:
lm(formula=sales~price)
Residuals:
Min
1QMedian
3Q
Max
**
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)(Intercept)<2e-16***price<2e-16***---Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
Residualstandarderror:
on102degreesoffreedom
MultipleR-squared:
0.7926,
AdjustedR-squared:
F-statistic:
on1and102DF,
p-value:
<
>acf(residuals(chip.m1),ci.type='ma')
由于回归后的残差自相关在四阶是显著的,因此我们要对其进行再一步的分析
**eacf(residuals(chip.m1))AR/MA
0xxxxooxxooooo1xooxooooooooo2xxoxooooooooo3xxoxooooooooo4oxxoooooooooo5xxxoxoooooooo6xxoxxxooooooo7xoxoooooooooo
oooooooo
Eacf
推荐其残差包含一个以〔
1,4〕为顶点为的零值三角形,从而说明其为
arma
〔1,4〕
模型,因此可将对数化销售量拟合成对于价格序列的带有
ARMA
〔1,4〕误差的回归模型。
模拟ARMA〔1,4〕初探
chip.m2=arima(sales,order=c(1,0,4),xreg=data.frame(price))
Call:
arima(x=sales,order=c(1,0,4),xreg=data.frame(price))
**
Coefficients:
ar1
ma1
ma2
ma3
ma4
intercept
price
s.e.
sigma^2estimatedas0.02556:
loglikelihood=42.35,
aic=
结果说明
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对应的系数。
>
Call:
arima(x=sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,0,NA,0,NA,NA,NA))
Coefficients:
ar1
ma10
ma2
ma30
ma4
intercept
price
**
s.e.
0
0
sigma^2estimatedas0.02572:
loglikelihood=42.09,aic=此模型的AR1系数项并不显著,所以再次调整模型
>chip.m4=arima(sales,order=c(0,0,4),xreg=data.frame(price),fixed=c(0,NA,0,NA,NA,NA))>
Call:
arima(x=sales,order=c(0,0,4),xreg=data.frame(price),fixed=c(0,NA,0,NA,NA,NA))
Coefficients:
ma1
ma2
ma3
ma4intercept
price
0
0
s.e.
0
0
sigma^2estimatedas0.02623:
loglikelihood=41.02,aic=此时模型建立完成,与一般线性回归比拟,两模型的截距项与价格项系数是相似的,但是用
时间序列估计的标准误差比用简单OLS回归所得的结果大约低10%,这说明了如下的结论,
**
即简单的OLS估计量具有一致性,但相关联的标准误差一般却是不可靠的。
对最终模型进行诊断分析
tsdiag(chip.m4)附m2=arima(days,order=c(0,0,2),xreg=data.frame(AO=seq(days)==129))#
拟合含有AO值时用xreg设置,假设无IO可直接用arima拟合m3=arimax(days,order=c(0,0,2),xreg=data.frame(AO=seq(days)==
129),io=c(63))#拟合含有IO值的要用arimax。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 时间 序列 回归 模型 实现