8时间序列回归模型R实现.docx
- 文档编号:7699604
- 上传时间:2023-01-25
- 格式:DOCX
- 页数:19
- 大小:341.65KB
8时间序列回归模型R实现.docx
《8时间序列回归模型R实现.docx》由会员分享,可在线阅读,更多相关《8时间序列回归模型R实现.docx(19页珍藏版)》请在冰豆网上搜索。
8时间序列回归模型R实现
时间序列回归模型
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')
>air.m1
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:
ma1sma1Dec96Jan97Dec02I911-MA0I911.1-AR1I911.1-MA0
-0.3825-0.64990.0989-0.06900.0810-0.09490.8139-0.2715
s.e.0.09260.11890.02280.02180.02020.04620.09780.0439
sigma^2estimatedas0.0006721:
loglikelihood=219.99,aic=-423.98
画图
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事件的影响效应才平息,航班客运量恢复了正常。
2异常值
在时间序列中异常有两种,可加异常和新息异常,分别记AO和IO。
2.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)
>y
TimeSeries:
Start=1
End=100
Frequency=1
[1]0.49180881-0.22323665-0.99151270-0.73387818-0.67750094-1.14472133-2.14844671-2.49530794
[9]-1.50355358-2.12615253-0.556517130.413263440.518691291.862106052.199354722.60210165
[17]0.791300030.262654262.934148573.990458893.608226781.17845765-0.87682948-1.20637799
[25]-1.39501221-0.188321711.229998271.468148502.666474913.234174692.603496241.49513215
[33]1.488521420.957392191.300116541.734440532.848251033.732146554.235794563.37049790
[41]2.027839551.41218929-0.29974176-1.58712591-1.340808780.107476091.446510811.67809487
[49]-0.34663129-0.502914590.01739605-0.014264740.942172040.39046221-0.398835301.60638918
[57]1.706682011.375181941.918245340.14254056-2.88169481-3.30372327-1.74068408-3.24868057
[65]-3.89415683-3.45920240-1.110420780.679597440.670510840.443940611.895360602.36063873
[73]2.005594430.864433240.468475720.723384981.602150981.259222771.531808590.96289779
[81]1.077121881.423863540.56318008-0.46689543-0.91861106-1.92947085-2.18188785-1.02759087
[89]2.310882723.138473193.012378813.434548072.315394942.449098732.915891411.12648908
[97]-0.081238710.444125790.26116418-0.45815484
>y[10]<-10
2.1.2模型初步判断
>acf(y)
>pacf(y)
>eacf(y)
AR/MA
012345678910111213
0xxoooooooooooo
1oooooooooooooo
2oooooooooooooo
3oxoooooooooooo
4oxoooooooooooo
5xxoooooooooooo
6xooooooooooooo
7oxoooooooooooo
从三个的结果来看,可以初步分析y是AR
(1)模型
2.1.3对模型时行拟合
>m1=arima(y,order=c(1,0,0))
>m1
Call:
arima(x=y,order=c(1,0,0))
Coefficients:
ar1intercept
0.54190.7096
s.e.0.08310.3603
2.1.4对模拟模型进行异常值探测
>detectAO(m1)
[,1][,2][,3]
ind9.00000010.00000011.000000
lambda2-4.0184129.068982-4.247367
>detectAO(m1,robust=F)
[,1]
ind10.000000
lambda27.321709
>detectIO(m1)
[,1][,2]
ind10.00000011.00000
lambda17.782013-4.67421
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))
>m2
Call:
arima(x=y,order=c(1,0,0),xreg=data.frame(AO=seq(y)==10))
Coefficients:
ar1interceptAO
0.80720.569810.9940
s.e.0.05700.51290.8012
sigma^2estimatedas1.059:
loglikelihood=-145.29,aic=296.58
>detectAO(m2)
[1]"NoAOdetected"
>detectIO(m2)
[1]"NoIOdetected"
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]"NoAOdetected"
>detectIO(m3)
[1]"NoIOdetected"
>tsdiag(m3)
>m3
Call:
arima(x=y,order=c(1,0,1),xreg=data.frame(AO=seq(y)==10))
Coefficients:
ar1ma1interceptAO
0.65960.61540.585011.1781
s.e.0.07990.07960.41320.4755
sigma^2estimatedas0.793:
loglikelihood=-131.16,aic=270.33
模型的拟合效果是显著提高。
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.co2
Call:
arima(x=co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
Coefficients:
ma1sma1
-0.5792-0.8206
s.e.0.07910.1137
sigma^2estimatedas0.5446:
loglikelihood=-139.54,aic=283.08
>detectAO(m1.co2)
[1]"NoAOdetected"
>detectIO(m1.co2)
[,1]
ind57.000000
lambda13.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:
ma1sma1IO-57
-0.5925-0.82742.6770
s.e.0.07750.10160.7246
sigma^2estimatedas0.4869:
loglikelihood=-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.2Log(销售量)与价格数据的相关性分析
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:
Min1QMedian3QMax
-0.54950-0.123730.006670.131360.45170
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)15.8900.21773.22<2e-16***
price-2.4890.126-19.75<2e-16***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.188on102degreesoffreedom
MultipleR-squared:
0.7926,AdjustedR-squared:
0.7906
F-statistic:
389.9on1and102DF,p-value:
<2.2e-16
>acf(residuals(chip.m1),ci.type='ma')
由于回归后的残差自相关在四阶是显著的,因此我们要对其进行再一步的分析
>eacf(residuals(chip.m1))
AR/MA
012345678910111213
0xxxxooxxoooooo
1xooxoooooooooo
2xxoxoooooooooo
3xxoxoooooooooo
4oxxooooooooooo
5xxxoxooooooooo
6xxoxxxoooooooo
7xoxooooooooooo
Eacf推荐其残差包含一个以(1,4)为顶点为的零值三角形,从而表明其为arma(1,4)模型,因此可将对数化销售量拟合成对于价格序列的带有ARMA(1,4)误差的回归模型。
4.2.3模拟ARMA(1,4)初探
>chip.m2=arima(sales,order=c(1,0,4),xreg=data.frame(price))
>chip.m2
Call:
arima(x=sales,order=c(1,0,4),xreg=data.frame(price))
Coefficients:
ar1ma1ma2ma3ma4interceptprice
0.1989-0.05540.25210.07350.526915.7792-2.4234
s.e.0.18430.16600.08650.10840.13760.21660.1247
sigma^2estimatedas0.02556:
loglikelihood=42.35,aic=-70.69
结果表明ma1,ma3的系数并不显著,即可认为其系数为0
4.2.4调整模型
>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.m3
Call:
arima(x=sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,
0,NA,0,NA,NA,NA))
Coefficients:
ar1ma1ma2ma3ma4interceptprice
0.144400.267600.521015.8396-2.4588
s.e.0.098500.085800.11710.20270.1166
sigma^2estimatedas0.02572:
loglikelihood=42.09,aic=-74.18
此模型的AR1系数项并不显著,所以再次调整模型
>chip.m4=arima(sales,order=c(0,0,4),xreg=data.frame(price),fixed=c(0,NA,0,NA,NA,NA))
>chip.m4
Call:
arima(x=sales,order=c(0,0,4),xr
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 时间 序列 回归 模型 实现