时间序列分析R语言程序.doc

1、1.5.csv,header=T)年中国年纱产量序列时序附录绘制1964 1999# 例2.1plot(Data1.5,type=o,xlim=c(1950,2010),ylim=c(60,100)）图（数据见附录1.2)Data1.2=read.csv(C:UsersAdministratorDesktoptdat1.5=Data1.5,2没有标题T 附录1.2.csv,header=T)# 如果有标题，用；a1.5=acf(tdat1.5) F用 plot(Data1.2,type=o)白噪声检验#Box.test(tdat1.5,type=Ljung-Box,lag=6)2.1 续 #例

2、Box.test(tdat1.5,type=Ljung-Box,lag=12) tdat1.2=Data1.2,2a1.2=acf(tdat1.2)模型拟合序列2.5 续选择合适的ARMA#例 acf(tdat1.5)pacf(tdat1.5)根据自相关系数图和偏自相关系数图可以判断为12 月平均每头奶#年例 #2.2绘制1962 年1 月至1975）模型牛产奶量序列时序图（数据见附录1.3）（1ARData1.3=read.csv(C:UsersAdministratorDesktop口径的求法在文档上P81 # 例 2.5续#P831.3.csv,header=F)附录然大似阵tdat1.

3、3=as.vector(t(as.matrix(Data1.3)1:168#矩转arima(tdat1.5,order=c(1,0,0),method=ML)#极 置转向量估计ar1=arima(tdat1.5,order=c(1,0,0),method=ML) plot(tdat1.3,type=l)summary(ar1)例 2.2 续 #ev=ar1$residualsacf(tdat1.3)#把字去掉acf(ev) pacf(tdat1.3)pacf(ev) 年北京市每年最高气温序19982.3# 例绘制 1949 列时序图Data1.4=read.csv(C:UsersAdminis

4、tratorDesktop #参数的显著性检验t1=0.6914/0.09891.4.csv,header=T)附录 p1=pt(t1,df=48,lower.tail=F)*2 plot(Data1.4,type=o)的显著性检验#ar1t2=81.5509/ 1.7453#不会定义坐标轴p2=pt(t2,df=48,lower.tail=F)*22.3# 例续tdat1.4=Data1.4,2残差白噪声检验#Box.test(ev,type=Ljung -Box,lag=6,fitdf=1) a1.4=acf(tdat1.4)Box.test(ev,type=Ljung -Box,lag=

5、12,fitdf=1)例 #2.3 续 Box.test(tdat1.4,type=Ljung-Box,lag=6)P94 续预测及置信区间 #例 2.5predict(arima(tdat1.5,order=c(1,0,0),n.ahead=5) Box.test(tdat1.4,type=Ljung -Box,lag=12) tdat1.5.fore=predict(arima(tdat1.5,order=c(1,0,0),n.ahead=5) 个服从标准正态分布的白噪声#随机产生2.41000 例序列观察值，并绘制时序图U=tdat1.5.fore$pred+1.96*tdat1.5.f

6、ore$se Data2.4=rnorm(1000,0,1)L=tdat1.5.fore$pred -1.96*tdat1.5.fore$se Data2.4plot(Data2.4,type=l)plot(c(tdat1.5,tdat1.5.fore$pred),type=l,col=1:2)#例续 2.4lines(U,col=lue,lty=dashed)a2.4=acf(Data2.4)lines(L,col=lue,lty=dashed)2.4 例 #续 Box.test(Data2.4,type=Ljung-Box,lag=6)Box.test(Data2.4,type=Ljung

7、-Box,lag=12)续例例3.53.53.1.1 # 例plot.ts(arima.sim(n=100,list(ar=0.8) #方法一#年北京市城乡居民定期储蓄1998 1950 对 2.5例方法二#x0=runif(1)所占比例序列的平稳性与纯随机性进行检验x=rep(0,1500)x1=x1x1=0.8*x0+rnorm(1)x2= -x1-0.5*x0+rnorm(1) for(i in 2:length(x)for(i in 3:length(x)xi=0.8*xi-1+rnorm(1)xi= -xi -1 -0.5*xi -2+rnorm(1) plot(x1:100,typ

8、e=l)plot(x1:100,type=l) acf(x)acf(x) pacf(x)pacf(x)#拟合图没有画出来均值和方差3.1.2 #例 x0=runif(1)smu=mean(x) x=rep(0,1500)svar=var(x) x1= -1.1*x0+rnorm(1)for(i in 2:length(x)xi= -1.1*xi -1+rnorm(1) 3.3例 AR （ 1）模型的方差plot(x1:100,type=l)acf(x)51 页 mvar=1/(1 -0.82) # 书上pacf(x)3.2#例求平稳mu=0总体均值方差) 3.1.3 # 例方法一 plot.t

9、s(arima.sim(n=100,list(ar=c(1,-0.5)样本均值方差#cat(sample mean)#方法二x0=runif(1)x1=runif(1) 3.4例题 #svar=(1+0.5)/(1 -0.5)*(1 -1-0.5)*(1+1 -0.5) x=rep(0,1500)x1=x1x2=x1 -0.5*x0+rnorm(1)自相关系数图截尾和偏自相关系模型例题 #3.6 MAfor(i in 3:length(x)数图拖尾 #3.6.1 xi=xi-1 -0.5*xi -2+rnorm(1)plot(x1:100,type=l)法一： x=arima.sim(n=10

10、00,list(ma=-2) acf(x)plot.ts(x,type=l) pacf(x)acf(x)pacf(x) 3.1.4 # 例 x0=runif(1)x1=runif(1)法二 x=rep(0:1000) x=rep(0,1500)for(i in 1:1000)x1=x1xi=rnormi-2*rnormi -1 x2=x1+0.5*x0+rnorm(1)plot(x,type=l) for(i in 3:length(x)acf(x) xi=xi-1+0.5*xi -2+rnorm(1)pacf(x) plot(x1:100,type=l)acf(x)#3.6.2 pacf(x

11、)法一：x=arima.sim(n=1000,list(ma=-0.5)又一个式子plot.ts(x,type=l) x0=runif(1)acf(x) x1=runif(1)pacf(x)x=rep(0,1500)天的 57ARMA模型拟合加油站#例 3.8 法二选择合适的x=rep(0:1000)序列OVERSHORTData1.6=read.csv(C:UsersAdministratorDesktopfor(i in 1:1000)1.6.csv,header=F)附录 tdat1.6=as.vector(t(as.matrix(Data1.6)1:57 xi=rnormi-0.5*r

12、normi -1plot(tdat1.6,type=o) plot(x,type=l)acf(x)acf(tdat1.6)pacf(x)pacf(tdat1.6) #closure# 错误于 rnormi :类别为的对象不可以取子把字去掉集 最小二乘arima(tdat1.6,order=c(0,0,1),method=CSS)#3.6.3估计ma1=arima(tdat1.6,order=c(0,0,1),method=CSS)法一： summary(ma1)x=arima.sim(n=1000,list(ma=c(-4/5,16/25)ev=ma1$residuals plot.ts(x,

13、type=l)acf(ev) acf(x)pacf(ev)pacf(x)= 1), method = c(0, 0, # 错误于法二 : arima(tdat1.6, order CSS) :x=rep(0:1000)#x 必需为数值for(i in 1:1000)年 1985 模型拟合 1880 #例 3.9 选择合适的ARMAxi=rnormi-4/5*rnormi -1+16/25*rnormi -2全球气温改变差值差分序列plot(x,type=l)没有数据 #acf(x)矩估计例 3.12#例 3.10例 3.11pacf(x)次化学反应的过程数xi = rnormi- 4/5 *

14、rnormi - 1 + 16/25 * #错误于 703.13 对等时间间隔的连续 #例 rnormi - 2 : 据进行拟合Data1.8=read.csv(C:UsersAdministratorDesktop #更换参数长度为零 1.8.csv,header=F)3.6#例续 64 根据书上页来判断附录tdat1.8=as.vector(t(as.matrix(Data1.8)1:70plot(tdat1.8,type=o)，（合 3.7#例拟 ARMA11型）模并直观观察该模型自相关，x(t) -0.5x(t -1)=u(t) -0.8*(u -1)）模型 3）例 3.16AR（ 系

15、数和偏自相关系数的拖尾性。 3.15AR例 3.14AR（2）例（ 3# #的预测法一：x0=runif(1)。 #如果考得话就先。 。 x=rep(0,1000)x1=0.5*x0+rnorm(1)-0.8*rnorm(1)4.1 线性拟合消费支出数据#例Data4.1=read.csv(C:UsersAdministratorDesktopfor(i in 2:length(x)xi=0.5*xi-1+rnorm(1) -0.8*rnorm(1) 4.1.csv,header=T)例题 tdat4.1=Data4.1,2 plot(x,type=l)plot(Data4.1,type=o)

16、 acf(x)t=1:40pacf(x)线性拟合 lm4.1=lm(tdat4.1t) # 图和书上不一样# summary(lm4.1) # 法二 #返回拟合参数的统计量 x=arima.sim(n=1000,list(ar=0.5,ma= -0.8)返回被估计的系数coef(lm4.1) #acf(x)返回模拟值 fit4.1=fitted(lm4.1) #pacf(x)返回残差值residuals(lm4.1) #画时序图plot(tdat4.1,type=o) #图和书上一样#m1=lm(yx,data=lindat) #一元线性回归拟合lines(fit4.1,col=ed) #画拟

17、合图summary(m1)#例4.2曲线拟合上海证劵交易所t=1:96Data1.9=read.csv(C:UsersAdministratorDesktopthat=1015.5222+20.9318*t1.9.csv,header=F)附录plot(1:length(tdat1.11),noseandat,p)转阵tdat1.9=as.vector(t(as.matrix(Data1.9)1:130#矩拟合图和原来的图画在一起lines(that,type=l)置转向量#plot(tdat1.9,type=l)残差检验#t=1:130计算残差ev=noseandat-that#evt2=t

18、2一道矩阵就出毛病残差图plot(ev) #m1.9=lm(tdat1.9t+t2) #t=97:108that=983.5601+21.5908*t简单移动平均法4.3# 例q=that*seaind x4.3=c(5,5.4,5.8,6.2)s=c(tdat1.11,q)x4.3plot(1:108,s,type=b) y4.3=filter(x4.3,rep(1/4,4),sides=1)abline(v=96) y4.3for(i in 1:3)x1=x1xi+1=0.25*xi+1+0.75*xi差分运算例5.1 #Data1.2=read.csv(C:UsersAdministra

19、torDesktopcolumns undefined 1)误于 #.data.frame(x, i + :错selected没有标题T ；用附录1.2.csv,header=T)# 如果有标题，F警告信息： :用 # 此外 x=Data1.2#In Ops.factor(left, right) : *对因子没有意义plot(x,type=o)dx=diff(x,2)指数平滑法4.4#例plot(dx,type=o)#做不出来北京市民车辆拥有量序列4.5 例 5.2 二阶差分#例#Data1.12=read.csv(C:UsersAdministratorDesktop略略# 1.12.cs

20、v,header=T)附录x=Data1.124.6#例季节效应分析plot(x,type=o)例 #4.7 中国社会消费品零售总额序列综合分析Data1.11=read.csv(C:UsersAdministratorDesktopT 1.11.csv,header=T) # 附录第一行是标签，所以是一届差分dx=diff(x,2)#plot(dx,type=b)二阶差分ddx=diff(x,2,lag=1,difference=2) #横向全部读取，纵向tdat1=as.matrix(Data1.11,2:9) #plot(ddx,type=l)92 读取至列tdat1.11=as.vec

21、tor(tdat1)5.2 又#画时序图，plot(1:length(tdat1.11),tdat1.11,type=o)#例Data1.12=read.csv(C:UsersAdministratorDesktop先是横坐标， 后是纵坐标1.12.csv,header=T)附录求总的均值md=mean(tdat1.11)#plot(Data1.12,2,type=l) mdaxis(1,at=c(1950,19999)seaind=apply(tdat1,1,mean)/md #求季节因子x=ts(Data1.12,2)seaindplot(seaind,type=b) #一阶差分季节指数图

22、#dx=diff(x,lag=1,differences=1)plot(Data1.12 -1,1,dx,type=o) 消除季节因子的影响 noseandat=tdat1.11/seaind # 二阶差分 消除季节因子之 plot(1:length(tdat1.11),noseandat,p) #d2x=diff(dx)#plot(Data1.12 -c(1,2),1,d2x,type=o)后的散点图lindat=data.frame(x=1:length(noseandat),y=noseandat)二阶差分 # d2x=diff(x,differences=2)Box.test(dx,l

23、ag=18,type=Ljung-Box) plot(Data1.12-c(1,2),1,d2x,type=o)pacf(dx)m1=arima(x,2,order=c(0,1,1),method=CSS)平均每头奶牛产奶量#例 5.3 跳步差分m2=arima(dx,order=c(0,1,1),method=CSS)Data1.13=read.csv(C:UsersAdministratorDesktopev1=m1$residuals 1.13.csv,header=F)附录ev2=m2$residuals tdat1=as.matrix(Data1.13)plot(ev1,type=l

24、) tdat=t(tdat1)plot(ev2,type=l) x=as.vector(tdat)acf(ev1) xpacf(ev1) plot(x,xaxt=n,type=o)acf(ev2) axis(1,at=seq(1,169,24),seq(1962,1976,2)pacf(ev2)dx=diff(x) # 一步差分plot(dx,type=o)检验残差的白噪Box.test(ev1,lag=5,type=Ljung-Box) #步差分声序列d12x=diff(dx,lag=12)#12Box.test(ev1,lag=11,type=Ljung-Box) plot(d12x,ty

25、pe=o)Box.test(ev1,lag=17,type=Ljung-Box)例5.4过差分#没做好做预测# 例 5.6续px=predict(m1,n.ahead=10)例5.5你和随机游走模型#plot(x,type=o,ylim=c(0,500)之间000r=rnorm(1000,sd=10)#以十位等差，在1： lines(x,1,x,2+1.96*sqrt(61.95)个数据随机抽取 100lines(x,1,x,2-1.96*sqrt(61.95)xt=cumsum(r) # 由随机游走公式的出的模型公式没画出来 14#图 5plot(xt,type=l)#随机游走的图形5.6

26、续 #dx=diff(xt) #做一阶差分例 m3=arima(x,2,order=c(0,1,1),method=ML)plot(dx,type=l) #一阶差分后的图形#m=mean(dx)均值p-171续 方差sd=var(dx) # #例5.6plot(x,xlim=c(1950,1990),ylim=c(0,300),type=o)用统计量检验随机#Box.test(dx,lag=12,type=Ljung)的函数,data=x) # 变量为时间tm1=lm(农业性年份？模型口径不会算summary(m1)# 自相关图acf(dx)#sj=arima(xt,order=c(0,1,0

27、) ,m1$fitted.value,col=red)年份lines(x$summary(sj)变量为一阶延迟#xt=x,2xy=xt -1又你和随机游走模型例#5.5xx=xt -length(xt) x=ts(cumsum(rnorm(1000,0,100)m2=lm(xyxx) ts.plot(x)summary(m2)m3=lm(xyxx+0)5.6 例 #对中国农业实际国民收入指数进行建模ARIMAsummary(m3)图份-1,m2$fitted.value,col=blue)#lines(x$ 年 29 5 1.14.csv,header=T)附录 x=Data1.14检验 #D

28、Wlibrary(lmtest)10 图 plot(x,type=o) #5dx=diff(x,2)dwtest(m2)# 加载程序包aa=dwtest(m1)plot(dx,type=o)Dh=(1 -aa$statistic/2)*sqrt(length(xt)-1)/(1 -(length) -1)*0.acf(dx)Box.test(dx,lag=6,type=Ljung-Box)009063 #Dh 统计量ev1=m1$residualsBox.test(dx,lag=12,type=Ljung-Box)dx1=diff(x)plot(ev1,type=o)plot(dx1,type

29、=o) m4=arima(ev1,order=c(2,0,0),fixed=c(NA,NA,0),trannsfdx=diff(dx1,lag=4)orm.pars=F)plot(dx,type=o,ylim=c(-2,2)ev2=m3$residualsBox.test(dx,lag=6,type=Ljung-Box) plot(ev2,type=o)Box.test(dx,lag=12,type=Ljung-Box)m5=arima(ev2,order=c(2,0,0),fixed=c(NA,0),trannsform.pars=F)差分 6 5 Box.test(dx,lag=18,ty

30、pe=Ljung-Box) # 表m6=arima(xt,order=c(0,1,1),xreg=1:length(xt),method=序列具有很强的相关信息CSS)#模型拟合 acf(dx) pacf(dx)m=arima(x,order=c(4,1,0),fixed=c(NA,0,0,NA),transform.pars=FALSE,include.mean=F,method=CSS)#模型拟合的不对coef(m)例 #5.7 ARIMA#不会 # 5.8 疏系数模型妇女参数估计与检验#例 ev1=m$residualsData1.15=read.csv(C:UsersAdministr

31、atorDesktopBox.test(ev,lag=4,type=Ljung-Box)1.15.csv,header=T)附录 Box.test(ev,lag=10,type=Ljung-Box) x=Data1.15plot(x,type=o)dx=diff(x,2)5.10 乘积季节模型#例Data1.17=read.csv(C:UsersAdministratorDesktopplot(dx,type=o)acf(dx) 1.17.csv,header=F)附录 tdat1.17=as.vector(t(as.matrix(Data1.17)1:408 pacf(dx)x=ts(tdat1.17,