应用时间序列分析第一次作业

1、应用时间序列分析第一次作业P131.2 附录B中的B1，B2分别是北京地区19852000年的月平均气温和降水量数据，其中缺少1989年的数据，B2还缺少1995年1月数据。（1） 用简单方法补齐1989年的数据和B2中1995年1月的数据，给出季节项的周期；（2） 对19902000年的两种数据个给出一种计算趋势项、季节项和随机项的公式；（3） 利用（2）的公式对所述的数据进行时间序列的分解计算，用数据图列出结果；（4） 用（2）中的结果补充1989年的数据解：（1）对于B1：使用分段趋势法有：（将趋势项定为年平均值）趋势项：B1减去趋势项剩下的部分为季节项和随机项。剩下的部分对每月取平均作

2、为季节项，则有：季节项：余下的部分则为随机项：对于B2：同理：（4）对于B1：1989年的数据为当年的趋势项加季节项对于B2：989年的数据为当月的趋势项加季节项1.3 附录B中的B6是19731978年美国在意外事故中的死亡人数。利用至少两种方法对该时间序列进行分解，要求如下：（1） 画出数据图，给出数据的周期T； 数据图如下： 使用R中decompose函数进行时序分解如下： 可以看出趋势项大致为二次曲线，因此可以考虑采用二次曲线来拟合趋势项。（2） 给出趋势项、季节项和随机项的计算公式； 方法一：二次曲线 认为趋势项满足二次线性方程，由最小二乘公式计算出趋势项如图： B6减去趋势项剩下的

3、部分为季节项和随机项。剩下的部分对每月取平均作为季节项，剩下的部分即为随机项。 方法二：回归方法（多元）认为趋势项满足多元线性方程，由最小二乘公式，计算出，再计算得出趋势项、季节项和随机项： （3） 画出趋势项、季节项和随机项的数据图； 方法一： 趋势项： 季节项： 随机项： 方法二： 趋势项、季节项： 随机项：（4） 对1979年的意外死亡人数做出预测。 方法一： 预测：1979年的数据为当月的趋势项加季节项： 预测图： 方法二： 预测：1979年的数据为回归方程结果： 预测图： 程序代码：#B1B1=matrix(c(-4.7,-3.7,-3.6,-2.9,-4.9,-2.3,-1.1,-

4、3.7,-1.6,-0.7,-2.2,-3.8,-3.9,-1.6,-6.4,-1.9,-1.8,0.1,-1.4,-0.6,0.1,1.8,1.6,0.8,2.1,-0.4,1.3,2.4,2.1,-1.5,3.4,6.9,4.1,4.4,7.6,4.4,6.7,8.1,5.6,7.7,6.2,8.7,7.6,4.7,8.0,14.8,15.0,13.5,15.0,13.7,13.9,15.5,14.0,17.3,14.7,14.3,14.5,15.0,14.4,14.6,19.5,21.3,19.9,20.1,19.6,19.9,20.5,21.5,21.0,19.8,21.6,20.0,

5、19.9,19.3,20.4,24.2,25.3,23.3,24.9,24.8,24.1,23.5,25.4,26.8,24.3,25.4,24.6,23.6,25.3,26.7,25.5,25.1,26.6,25.8,25.6,25.9,26.8,25.2,27.7,25.9,25.5,28.2,26.5,28.0,29.6,25.0,24.5,24.8,24.4,25.4,27.1,24.6,25.2,26.5,25.4,23.9,26.6,25.1,25.5,25.7,18.6,19.8,21.0,21.2,20.2,20.4,20.5,21.3,21.1,19.0,20.7,18.6,

6、22.2,20.9,21.8,13.8,11.4,13.7,14.1,15.3,13.8,12.2,13.9,14.1,14.5,12.8,14.0,14.8,12.9,12.6,3.8,3.4,3.9,6.9,6.4,4.6,3.4,3.7,6.4,7.7,4.2,5.4,4.0,5.9,3.0,-3.6,-1.7,-0.3,-0.2,-0.8,-1.8,-0.3,-0.8,-1.4,-0.4,0.9,-1.5,0.1,-0.7,-0.6),ncol=12)B1.mean=apply(B1,2,mean)B1=rbind(B11:4,B1.mean,B15:15,deparse.level

7、=0)B1.trend=apply(B1,1,mean)#trendplot(rep(c(t(B1.trend),rep(12,16),type="l",ylab="trend")B1.temp=B1-B1.trendB1.season=apply(B1.temp,2,mean)B1.season=B1.season-sum(B1.season)/12#seasonalplot(rep(c(t(B1.season),16),type="l",ylab="seasonal")B1.rand=t(t(B1.temp)-

8、B1.season)#randomplot(c(t(B1.rand),type="h",ylab="random")# transform matrix into time series B1=c(t(B1)B1.ts=ts(data=B1,start=1985,frequency=12)plot(B1.ts)# use decompose functionDecom.B1=decompose(B1.ts)plot(Decom.B1)#B2B2=c(1.5,7.5,7.8,13.6,24.5,32.0,289.5,297.7,38.4,3.8,4.6,0

9、.1,0.0,6.0,17.0,1.0,5.0,203.0,163.0,143.0,114.0,4.0,6.0,4.0,4.3,2.4,13.0,41.8,64.6,91.2,130.9,246.5,46.2,4.1,35.4,3.5,0.9,1.3,8.9,8.2,37.4,61.8,278.7,204.0,48.8,22.8,0.0,0.5,4.7,21.6,40.5,59.7,119.6,4.0,223.0,157.0,63.1,0.3,3.6,0.2,0.3,0.8,25.1,17.1,214.6,236.3,198.0,124.7,72.0,12.2,1.0,4.7,0.7,0.0,

10、3.4,10.5,52.8,69.4,153.9,141.4,54.5,38.1,16.7,0.1,3.7,1.5,0.3,16.9,8.6,39.2,206.4,158.5,18.3,9.9,43.4,0.0,0.0,5.0,0.0,1.9,66.0,23.6,459.2,214.2,15.2,10.3,12.7,5.1,NA,1.7,6.6,5.3,45.6,68.9,195.6,119.9,116.3,9.6,0.2,2.8,0.2,0.0,11.0,6.2,1.8,55.1,307.4,250.0,32.9,30.8,2.6,2.9,4.9,0.0,10.6,17.4,41.5,35.

11、5,139.8,83.2,44.1,43.0,2.1,8.8,1.3,26.3,4.3,54.7,61.5,142.9,247.9,114.4,4.7,61.8,11.3,0.6,0.0,0.0,5.2,33.6,32.4,23.8,62.7,63.5,44.5,3.9,9.5,0.7,11.9,0.0,8.8,18.3,37.7,19.0,61.5,150.5,18.4,35.2,9.7,0.1)B2=matrix(B2,ncol=12,byrow=TRUE)#fill the missing valueB2.mean=apply(B2,2,mean,na.rm=TRUE)B2=rbind(

12、B21:4,B2.mean,B25:15,)B2=c(t(B2)B2is.na(B2)=B2.mean1# transform matrix into time series B2.ts=ts(data=B2,start=1985,frequency=12)plot(B2.ts)# use decompose functionDecom.B2=decompose(B2.ts)plot(Decom.B2)# least-squares fittingY=matrix(c(rep(1,192),1:192),ncol=192,byrow=TRUE)Result=solve(Y%*%t(Y)%*%Y

13、%*%B2B2.trend=c(t(Y)%*%Result)#trendplot(B2.trend,type="l",ylab="trend")B2.temp=B2-B2.trendB2.temp=matrix(B2.temp,ncol=12,byrow=TRUE)B2.season=apply(B2.temp,2,mean)B2.season=B2.season-sum(B2.season)/12#seasonalplot(rep(c(t(B2.season),16),type="l",ylab="seasonal&quo

14、t;)B2.rand=t(t(B2.temp)-B2.season)#randomplot(c(t(B2.rand),type="h",ylab="random")#(4)#B1B1.1989=B1.trend5+B1.season#B2B2.1989=B2.trend(4*12+1):(4*12+12)+B2.season#Question 1.3B6=c(9007,8106,8928,9137,10017,10826,11317,10744,9713,9938,9161,8927,7750,6981,8038,8422,8714,9512,10120

15、,9823,8743,9192,8710,8680,8162,7306,8124,7870,9387,9556,10093,9620,8285,8433,8160,8034,7717,7461,7776,7925,8634,8945,10078,9179,8037,8488,7874,8647,7792,6957,7726,8106,8890,9299,10625,9302,8314,8850,8265,8796,7836,6892,7791,8129,9115,9434,10484,9827,9110,9070,8633,9240)B6=matrix(B6,ncol=12,byrow=TRU

16、E)B6.mean=apply(B6,2,mean)B6=c(t(B6)# transform matrix into time series B6.ts=ts(data=B6,start=1973,frequency=12)plot(B6.ts)# use decompose functionDecom.B6=decompose(B6.ts)plot(Decom.B6)# least-squares fitting(quadratic)Y=matrix(c(rep(1,72),1:72,1:72*1:72),ncol=72,byrow=TRUE)Result=solve(Y%*%t(Y)%*

17、%Y%*%B6B6.trend=c(t(Y)%*%Result)#trendplot(B6.trend,type="l",ylab="trend")B6.temp=B6-B6.trendB6.temp=matrix(B6.temp,ncol=12,byrow=TRUE)B6.season=apply(B6.temp,2,mean)B6.season=B6.season-sum(B6.season)/12#seasonalplot(rep(c(t(B6.season),6),type="l",ylab="seasonal&qu

18、ot;)B6.rand=t(t(B6.temp)-B6.season)#randomplot(c(t(B6.rand),type="h",ylab="random")# predict 1979Y.1979=matrix(c(rep(1,12),73:84,73:84*73:84),ncol=12,byrow=TRUE)B6.1979_1=c(t(Y.1979)%*%Result)+B6.seasonplot(1:72,B6,type="l",xlim=c(0,84),main="Prediction",xlab="Time")lines(73:84,B6.1979_1,col="BLUE",lty=3)# least-squ