應(yīng)用時(shí)間序列分析第一次作業(yè).doc

應(yīng)用時(shí)間序列分析第一次作業(yè)P131.2 附錄B中的B1，B2分別是北京地區(qū)19852000年的月平均氣溫和降水量數(shù)據(jù)，其中缺少1989年的數(shù)據(jù)，B2還缺少1995年1月數(shù)據(jù)。（1） 用簡(jiǎn)單方法補(bǔ)齊1989年的數(shù)據(jù)和B2中1995年1月的數(shù)據(jù)，給出季節(jié)項(xiàng)的周期；（2） 對(duì)19902000年的兩種數(shù)據(jù)個(gè)給出一種計(jì)算趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和隨機(jī)項(xiàng)的公式；（3） 利用（2）的公式對(duì)所述的數(shù)據(jù)進(jìn)行時(shí)間序列的分解計(jì)算，用數(shù)據(jù)圖列出結(jié)果；（4） 用（2）中的結(jié)果補(bǔ)充1989年的數(shù)據(jù)解：（1）對(duì)于B1：使用分段趨勢(shì)法有：（將趨勢(shì)項(xiàng)定為年平均值）趨勢(shì)項(xiàng)：B1減去趨勢(shì)項(xiàng)剩下的部分為季節(jié)項(xiàng)和隨機(jī)項(xiàng)。剩下的部分對(duì)每月取平均作為季節(jié)項(xiàng)，則有：季節(jié)項(xiàng)：余下的部分則為隨機(jī)項(xiàng)：對(duì)于B2：同理：（4）對(duì)于B1：1989年的數(shù)據(jù)為當(dāng)年的趨勢(shì)項(xiàng)加季節(jié)項(xiàng)對(duì)于B2：989年的數(shù)據(jù)為當(dāng)月的趨勢(shì)項(xiàng)加季節(jié)項(xiàng)1.3 附錄B中的B6是19731978年美國(guó)在意外事故中的死亡人數(shù)。利用至少兩種方法對(duì)該時(shí)間序列進(jìn)行分解，要求如下：（1） 畫出數(shù)據(jù)圖，給出數(shù)據(jù)的周期T； 數(shù)據(jù)圖如下： 使用R中decompose函數(shù)進(jìn)行時(shí)序分解如下： 可以看出趨勢(shì)項(xiàng)大致為二次曲線，因此可以考慮采用二次曲線來(lái)擬合趨勢(shì)項(xiàng)。（2） 給出趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和隨機(jī)項(xiàng)的計(jì)算公式； 方法一：二次曲線 認(rèn)為趨勢(shì)項(xiàng)滿足二次線性方程，由最小二乘公式計(jì)算出趨勢(shì)項(xiàng)如圖： B6減去趨勢(shì)項(xiàng)剩下的部分為季節(jié)項(xiàng)和隨機(jī)項(xiàng)。剩下的部分對(duì)每月取平均作為季節(jié)項(xiàng)，剩下的部分即為隨機(jī)項(xiàng)。 方法二：回歸方法（多元）認(rèn)為趨勢(shì)項(xiàng)滿足多元線性方程，由最小二乘公式，計(jì)算出，再計(jì)算得出趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和隨機(jī)項(xiàng)： （3） 畫出趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和隨機(jī)項(xiàng)的數(shù)據(jù)圖； 方法一： 趨勢(shì)項(xiàng)： 季節(jié)項(xiàng)： 隨機(jī)項(xiàng)： 方法二： 趨勢(shì)項(xiàng)、季節(jié)項(xiàng)： 隨機(jī)項(xiàng)：（4） 對(duì)1979年的意外死亡人數(shù)做出預(yù)測(cè)。 方法一： 預(yù)測(cè)：1979年的數(shù)據(jù)為當(dāng)月的趨勢(shì)項(xiàng)加季節(jié)項(xiàng)： 預(yù)測(cè)圖： 方法二： 預(yù)測(cè)：1979年的數(shù)據(jù)為回歸方程結(jié)果： 預(yù)測(cè)圖： 程序代碼：#B1B1=matrix(c(-4.7,-3.7,-3.6,-2.9,-4.9,-2.3,-1.1,-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,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,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 =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)-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.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,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.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(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%*%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)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,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=TRUE)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)%*%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)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-squares fitting(multivariate)Y2=matrix(rep(diag(12),6),ncol=72)Y2=rbind(t(1:72),Y2,deparse.level =0)Result2=solve(Y2%*%t(Y2)%*%Y2%*%B6B6.trend2=c(t(Y2)%*%Result2)#trend+seasonalplot(B6.trend2,type=l,ylab=trend2)B6.rand2=B6-B6.tre