时间序列分析R语言程序

1、#例2.1绘制1964 1999年中国年纱产量序列时序 图(数据见附录1.2)Data1.2=read.csv(C:UsersAdmi nistratorWDesktop附录1.2.csv,header=T)#如果有标题，用T;没有标 题用Fplot(Data1.2,type=o)#例2.1续tdat1.2=Data1.2,2 a1.2=acf(tdat1.2)#例2.2绘制1962年1月至1975年12月平均每头奶牛产奶量序列时序图(数据见附录1.3)Data1.3=read.csv(C:UsersAdmi nistratorWDesktop附录 1.3.csv,header=F) tdat

2、1.3=as.vector(t(as.matrix(Data1.3)1:168# 矩 阵 转置转向量plot(tdat1.3,type=T)#例2.2续acf(tdat1.3)# 把字去掉pacf(tdat1.3)#例2.3绘制1949 1998年市每年最高气温序列时 序图Data1.4=read.csv(C:UsersAdmi nistratorWDesktop附录 1.4.csv,header=T)plot(Data1.4,type=o)#不会定义坐标轴#例2.3续tdat1.4=Data1.4,2a1.4=acf(tdat1.4)#例2.3续Box.test(tdat1.4,type=L

3、ju ng-Box,lag=6)Box.test(tdat1.4,type=Lju ng-Box,lag=12)#例2.4随机产生1000个服从标准正态分布的白噪声 序列观察值，并绘制时序图Data2.4=rnorm(1000,0,1)Data2.4plot(Data2.4,type=T)#例2.4续a2.4=acf(Data2.4)#例2.4续Box.test(Data2.4,type=Lju ng-Box,lag=6)Box.test(Data2.4,type=Lju ng-Box,lag=12)#例2.5对1950 1998年市城乡居民定期储蓄所占 比例序列的平稳性与纯随机性进行检验Da

4、ta1.5=read.csv(C:UsersAdmi nistratorWDesktop附录 1.5.csv,header=T)plot(Data1.5,type=o,xlim=c(1950,2010),ylim=c(60,100)tdat1.5=Data1.5,2a1.5=acf(tdat1.5)#白噪声检验Box.test(tdat1.5,type=Lju ng-Box,lag=6)Box.test(tdat1.5,type=Lju ng-Box,lag=12)#例2.5续选择合适的 ARMA模型拟合序列acf(tdat1.5)pacf(tdat1.5)#根据自相关系数图和偏自相关系数图可

5、以判断为AR(1)模型#例2.5续P81 口径的求法在文档上#P83arima(tdat1.5,order=c(1,0,0),method=ML)# 极大似 然估计ar1= arima(tdat1.5,order=c(1,0,0),method=ML) summary(ar1)ev=ar1$residualsacf(ev)pacf(ev)#参数的显著性检验t1=0.6914/0.0989p1=pt(t1,df=48,lower.tail=F)*2#ar1的显著性检验t2=81.5509/ 1.7453p2=pt(t2,df=48,lower.tail=F)*2#残差白噪声检验Box.test(

6、ev,type=Lj un g-Box,lag=6,fitdf=1)Box.test(ev,type=Lj un g-Box,lag=12,fitdf=1)#例2.5续P94预测及置信区间predict(arima(tdat1.5,order=c(1,0,0), n. ahead=5)tdat1.5.fore=predict(arima(tdat1.5,order=c(1,0,0), n.a head=5)U=tdat1.5.fore$pred+1.96*tdat1.5.fore$seL=tdat1.5.fore$pred-1.96*tdat1.5.fore$seplot(c(tdat1.5,

7、tdat1.5.fore$pred),type=T,col=1:2)lin es(U,col=blue,lty=dashed)lin es(L,col=blue,lty=dashed)#例3.1.1例3.5 例3.5续#方法一 plot.ts(arima.sim(n=100,list(ar=0.8)#方法二xO=ru nif(1)x=rep(0,1500)x1=0.8*x0+rnorm(1)for(i in 2:le ngth(x)xi=0.8*xi-1+rnorm(1) plot(x1:100,type=T) acf(x)pacf(x)#拟合图没有画出来#例 3.1.2xO=ru nif(1

8、)x=rep(0,1500)x1=-1.1*xO+rnorm(1)for(i in 2:le ngth(x)xi=-1.1*xi-1+rnorm(1) plot(x1:100,type=T) acf(x)pacf(x)#例 3.1.3方法一plot.ts(arima.sim( n=100,list(ar=c(1,-0.5)#方法二xO=ru nif(1)x1=ru nif(1)x=rep(0,1500)x1=x1x2=x1-0.5*x0+rnorm(1)for(i in 3:le ngth(x)xi=xi-1-0.5*xi-2+rnorm(1) plot(x1:100,type=T) acf(

9、x)pacf(x)#例 3.1.4xO=ru nif(1)x1=ru nif(1)x=rep(0,1500)x1=x1x2=x1+0.5*x0+rnorm(1)for(i in 3:le ngth(x)xi=xi-1+0.5*xi-2+rnorm(1)plot(x1:100,type=T)acf(x)pacf(x)又一个式子xO=ru nif(1)x1=ru nif(1)x=rep(0,1500)x1=x1x2=-x1-0.5*x0+rnorm(1)for(i in 3:le ngth(x)xi=-xi-1-0.5*xi-2+rnorm(1)plot(x1:100,type=T)acf(x)p

10、acf(x)#均值和方差smu=mea n(x)svar=var(x)#例3.2求平稳AR( 1)模型的方差例3.3mu=0mvar=1/(1-0.8A2) # 书上 51 页#总体均值方差cat(populati on mean and var are,c(mu,mvar),n)#样本均值方差cat(sample mean and var are,c(mu,mvar),n)#例题3.4svar=(1+0.5)/(1-0.5)*(1-1-0.5)*(1+1-0.5)#例题3.6 MA模型 自相关系数图截尾和偏自相关系 数图拖尾#3.6.1法一：x=arima.sim (n=1000,list(

11、ma=-2)plot.ts(x,type=T)acf(x)pacf(x)法二x=rep(0:1000)for(i in 1:1000)xi=rnormi-2*rnormi-1 plot(x,type=T)acf(x)pacf(x) #3.6.2法一：x=arima.sim (n=1000,list(ma=-0.5)plot.ts(x,type=T)acf(x)pacf(x)法二x=rep(0:1000)for(i in 1:1000)xi=rnormi-0.5*rnormi-1plot(x,type=T)acf(x)pacf(x)#错误于rnormi:类别为closure的对象不可以取子集#3

12、.6.3法一：x=arima.sim( n=1000,list(ma=c(-4/5,16/25)plot.ts(x,type=T)acf(x)pacf(x)法二：x=rep(0:1000)for(i in 1:1000)xi=rnormi-4/5*rnormi-1+16/25*rnormi-2plot(x,type=T)acf(x)pacf(x)# 错误于 xi = rnormi - 4/5 * rnormi - 1 + 16/25 *rno rmi - 2:#更换参数长度为零#例3.6续根据书上64页来判断#例 3.7 拟合 ARMA ( 1，1) 模型，x(t)-0.5x(t-1)=u(t

13、)-0.8*(u-1)，并直观观察该模型自相 关系数和偏自相关系数的拖尾性。#法一：x0=ru nif(1)x=rep(0,1000)x1=0.5*x0+rnorm(1)-0.8*rnorm(1)for(i in 2:le ngth(x)xi=0.5*xi-1+rnorm(1)-0.8*rnorm(1) plot(x,type=T)acf(x)pacf(x)#图和书上不一样#法二x=arima.sim( n=1000,list(ar=0.5,ma=-0.8)acf(x)pacf(x)#图和书上一样#例3.8选择合适的 ARMA模型拟合加油站57天的OVERSHORT 序列Data1.6=rea

14、d.csv(C:UsersAdmi nistratorWDesktop 附录 1.6.csv,header=F)tdat1.6=as.vector(t(as.matrix(Data1.6)1:57 plot(tdat1.6,type=o)acf(tdat1.6)pacf(tdat1.6) # 把字去掉arima(tdat1.6,order=c(0,0,1),method=CSS)# 最小 二 乘估计ma1=arima(tdat1.6,order=c(0,0,1),method=CSS) summary(ma1)ev=ma1$residualsacf(ev)pacf(ev)# 错误于 arima

15、(tdat1.6, order = c(0, 0, 1), method =CSS):#x必需为数值#例3.9选择合适的 ARMA模型拟合1880 1985年全球气温改变差值差分序列#没有数据#例3.10例3.11例3.12#矩估计#例3.13对等时间间隔的连续70次化学反应的过程数据进行拟合Data1.8=read.csv(C:UsersAdmi nistratorWDesktop附录 1.8.csv,header=F)tdat1.8=as.vector(t(as.matrix(Data1.8)1:70plot(tdat1.8,type=o)#例 3.14AR (2 )例 3.15AR (3

16、 )例 3.16AR ( 3)模型 的预测#如果考得话就先。#例4.1线性拟合消费支出数据Data4.仁read.csv(C:UsersAdmi nistratorWDesktop例题 4.1.csv,header=T)tdat4.1=Data4.1,2 plot(Data4.1,type=o) t=1:40lm4.仁 Im(tdat4.1t) # 线性拟合 summary(lm4.1) #返回拟合参数的统计量 coef(lm4.1) #返回被估计的系数 fit4.1=fitted(lm4.1) # 返回模拟值 residuals(lm4.1) # 返回残差值 plot(tdat4.1,typ

17、e=o) #画时序图lines(fit4.1,col=red) # 画拟合图#例4.2曲线拟合证券交易所Data1.9=read.csv(C:UsersAdmi nistratorWDesktop附录 1.9.csv,header=F) tdat1.9=as.vector(t(as.matrix(Data1.9)1:130# 矩 阵 转置转向量plot(tdat1.9,type=T)t=1:130t2=tA2m1.9=lm(tdat1.9t+t2) #一道矩阵就出毛病#例4.3简单移动平均法x4.3=c(5,5.4,5.8,6.2)x4.3y4.3=filter(x4.3,rep(1/4,4)

18、,sides=1)y4.3for(i in 1:3)x1=x1xi+1=0.25*xi+1+0.75*xi# 错误于、.data.frame(x, i + 1) : undefined columns selected#此外：警告信息：#In Ops.factor(left, right) : * 对因子没有意义#例4.4指数平滑法#做不出来# 例 4.5#略略#例4.6季节效应分析#例4.7综合分析中国社会消费品零售总额序列Datal.1 仁read.csv(C:UsersAdmi nistratorWDeskto附录1.11.csv,header=T) #第一行是标签，所以是Ttdat仁a

19、s.matrix(Data1.11,2:9)# 横向全部读取，纵向读取2至9列tdat1.11=as.vector(tdat1)plot(1:length(tdat1.11),tdat1.11,type=o)# 画时序图，先是横坐标，后是纵坐标md=mea n(tdat1.11)# 求总的均值mdseai nd=apply(tdat1,1,mea n)/md #求季节因子sea indplot(seai nd,type=b) #季节指数图nosea ndat=tdat1.11/seai nd #消除季节因子的影响plot(1:le ngth(tdatl.ll), nosea ndat,p) #

20、消除季节因子之后的散点图lin dat=data.frame(x=1:le ngth( no sea ndat),y=no sea ndat)m仁lm(yx,data=li ndat) #一元线性回归拟合summary(ml)t=1:96that=1015.5222+20.9318*tplot(1:le ngth(tdatl.ll), nosea ndat,p)lin es(that,type=T)#拟合图和原来的图画在一起#残差检验ev=no sea ndat-that# 计算残差evplot(ev) #残差图t=97:108that=983.5601+21.5908*tq=that*sea

21、 ind s=c(tdat1.11,q)plot(1:108,s,type=b)abli ne(v=96)#例5.1差分运算Data1.2=read.csv(C:UsersAdmi nistratorWDesktop附录1.2.csv,header=T)#如果有标题，用T;没有标 题用Fx=Data1.2plot(x,type=o)dx=diff(x,2)plot(dx,type=o)#例5.2二阶差分市民车辆拥有量序列Data1.12=read.csv(C:UsersAdmi nistratorWDeskto 附录 1.12.csv,header=T)x=Data1.12plot(x,typ

22、e=o)dx=diff(x,2)#一届差分plot(dx,type=b)ddx=diff(x,2,lag=1,differe nce=2) #二阶差分plot(ddx,type=T)#例5.2又Data1.12=read.csv(C:UsersAdmi nistratorWDeskto 附录 1.12.csv,header=T)plot(Data1.12,2,type=T)axis(1,at=c(1950,19999)x=ts(Data1.12,2)dx=diff(x,lag=1,differe nces=1) # 一阶差分 plot(Data1.12-1,1,dx,type=o)d2x=di

23、ff(dx)# 二阶差分plot(Data1.12-c(1,2),1,d2x,type=o) d2x=diff(x,differe nces=2)# 二阶差分plot(Data1.12-c(1,2),1,d2x,type=o)#例5.3跳步差分平均每头奶牛产奶量Data1.13=read.csv(C:UsersAdmi nistratorWDeskto 附录 1.13.csv,header=F)tdat 仁as.matrix(Data1.13)tdat=t(tdat1)x=as.vector(tdat)xplot(x,xaxt= n,type=o) axis(1,at=seq(1,169,24

24、),seq(1962,1976,2) dx=diff(x) # 一步差分plot(dx,type=o)d12x=diff(dx,lag=12) #12 步差分 plot(d12x,type=o)#例5.4过差分#例5.5你和随机游走模型r=rnorm(1000,sd=10)# 以十位等差，在 1 : 000 之间随机抽取100个数据xt=cumsum(r) #由随机游走公式的出的模型公式plot(xt,type=l)# 随机游走的图形dx=diff(xt) # 做一阶差分plot(dx,type=T) # 一阶差分后的图形 m=mea n( dx)# 均值sd=var(dx) # 方差Box.

25、test(dx,lag=12,type=Lj un g)# 用统计量检验随机性acf(dx) #自相关图sj=arima(xt,order=c(0,1,0)summary(sj)#例5.5你和随机游走模型又x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)#例5.6对中国农业实际国民收入指数进行建模ARIMAData1.14=read.csv(C:UsersAdmi nistratorWDeskto 附录 1.14.csv,header=T)x=Data1.14plot(x,type=o) # 图 510dx=diff(x,2)plot(dx,type=o)acf(

26、dx)Box.test(dx,lag=6,type=Lj un g-Box)Box.test(dx,lag=12,type=Lju ng-Box)Box.test(dx,lag=18,type=Lju ng-Box)pacf(dx)m仁 arima(x,2,order=c(0,1,1),method=CSS) m2=arima(dx,order=c(0,1,1),method=CSS) ev1=m1$residualsev2=m2$residualsplot(ev1,type=T)plot(ev2,type=T)acf(ev1)pacf(ev1)acf(ev2)pacf(ev2)Box.tes

27、t(ev1,lag=5,type=Lju ng-Box)# 检验残差的白噪声序列Box.test(ev1,lag=11,type=Lju ng-Box)Box.test(ev1,lag=17,type=Lju ng-Box)#例5.6续 做预测#没做好px=predict(m1, n. ahead=10) plot(x,type=o,ylim=c(0,500) lin es(x,1,x,2+1.96*sqrt(61.95) lin es(x,1,x,2-1.96*sqrt(61.95)#图514没画出来#例5.6续 m3=arima(x,2,order=c(0,1,1),method=ML)#

28、例 5.6 续 p-171 plot(x,xlim=c(1950,1990),ylim=c(0,300),type=o) m1=lm(农业年份,data=x) #变量为时间t的函数 summary(m1) # ?模型口径不会算 lin es(x$ 年份，m1$fitted.value,col=red)#变量为一阶延迟xt=x,2xy=xt-1xx=xt-le ngth(xt)m2=lm(xyxx)summary(m2) m3=lm(xyxx+0)summary(m3)lines(x$ 年份-1,m2$fitted.value,col=blue) # 图529#DW检验library(lmtes

29、t)dwtest(m2)#加载程序包aa=dwtest(m1)Dh=(1-aa$statistic/2)*sqrt(le ngth(xt)-1)/(1-(le ngth)-1 )*0.009063 #Dh 统计量ev1=m1$residualsplot(ev1,type=o) m4=arima(ev1,order=c(2,0,0),fixed=c(NA,NA,0),tra nn sform.pars=F)ev2=m3$residualsplot(ev2,type=o) m5=arima(ev2,order=c(2,0,0),fixed=c(NA,0),tra nn sfor m.pars=F)

30、m6=arima(xt,order=c(0,1,1),xreg=1:le ngth(xt),method =CSS)#例 5.7 ARIMA#例5.8疏系数模型妇女Data1.15=read.csv(C:UsersAdmi nistratorWDeskto附录 1.15.csv,header=T)x=Data1.15plot(x,type=o)dx=diff(x,2)plot(dx,type=o)acf(dx)pacf(dx) m=arima(x,2,order=c(4,1,0),fixed=c(NA,0,0,NA),tra ns form.pars=FALSE,method=ML)# 不知道

31、疏系数模型是怎么判断的summary(m)ev=m$residualsBox.test(ev,lag=6,type=Ljung-Box,fitdf=2) # 阶数为原来的阶数减去参数的个数Box.test(ev,lag=12,type=Lj un g-Box,fitdf=2)Box.test(ev,lag=18,type=Ljung-Box,fitdf=2)# 结果和书上不一样#参数显著性检验t1=0.2583/0.11592.228645t2=0.3408/0.12252.782041 #t 统计量p仁pt(2.228645,df=57,lower.tail=F)*2 p2=pt(2.782041,df=57,lower.tail=F)*2 #p 值不大对#例5.9简单季节模型德国工人失业率Data1.16=read.csv(C:UsersAdmi nistratorWDeskto 附录 1.16.csv,header=F)tdat1.16=as.vector(t(as.matrix(Data1.16)1:120#绘制时序图x=ts(tdat1.1