8时间序列回归模型R实现Word格式.docx
- 文档编号:8184929
- 上传时间:2023-05-10
- 格式:DOCX
- 页数:20
- 大小:343.91KB
8时间序列回归模型R实现Word格式.docx
《8时间序列回归模型R实现Word格式.docx》由会员分享,可在线阅读,更多相关《8时间序列回归模型R实现Word格式.docx(20页珍藏版)》请在冰点文库上搜索。
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))
ar1intercept
0.54190。
7096
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。
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
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。
〉detectAO(m3)
[1]”NoAOdetected”
〉detectIO(m3)
[1]"
NoIOdetected”
〉tsdiag(m3)
〉m3
arima(x=y,order=c(1,0,1),xreg=data。
ar1ma1interceptAO
0.65960.61540.585011。
1781
e.0。
07990.07960。
41320。
4755
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
arima(x=co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
ma1sma1
5792—0。
8206
s.e.0。
07910.1137
sigma^2estimatedas0.5446:
loglikelihood=—139.54,aic=283。
08
〉detectAO(m1。
co2)
[1]”NoAOdetected”
〉detectIO(m1。
[,1]
ind57。
lambda13.752715
拟合含有新息异常的模型
〉m4。
co2=arimax(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),io=c(57))
arimax(x=co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),
io=c(57))
ma1sma1IO—57
-0。
5925—0。
82742。
6770
07750.10160。
7246
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)
lm(formula=sales~price)
Residuals:
Min1QMedian3QMax
—0.54950-0.123730。
006670.131360.45170
EstimateStd。
ErrortvaluePr(〉|t|)
(Intercept)15。
8900。
21773.22〈2e—16***
price-2.4890.126-19。
75〈2e-16***
-——
Signif。
codes:
0‘***'
001‘**'
01‘*’0.05‘.’0。
1‘’1
Residualstandarderror:
0.188on102degreesoffreedom
MultipleR—squared:
7926,AdjustedR—squared:
7906
F—statistic:
389.9on1and102DF,p—value:
〈2.2e-16
〉acf(residuals(chip。
m1),ci。
type='
ma’)
由于回归后的残差自相关在四阶是显著的,因此我们要对其进行再一步的分析
〉eacf(residuals(chip。
m1))
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))
m2
arima(x=sales,order=c(1,0,4),xreg=data.frame(price))
ar1ma1ma2ma3ma4interceptprice
1989—0。
05540.25210。
07350。
526915。
7792-2.4234
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
arima(x=sales,order=c(1,0,4),xreg=data。
frame(price),fixed=c(NA,
0,NA,0,NA,NA,NA))
0.144400。
267600.521015.8396-2.4588
098500。
085800。
11710。
20270。
1166
sigma^2estimatedas0.02572:
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))
m4
arima(x=sales,order=c(0,0,4),xreg=d
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 时间 序列 回归 模型 实现
![提示](https://static.bingdoc.com/images/bang_tan.gif)