時間序列分析R語言程序.doc

1、1.5.csv,header=T)年中國年紗產(chǎn)量序列時序附錄繪制1964 1999# 例2.1plot(Data1.5,type=o,xlim=c(1950,2010),ylim=c(60,100)）圖（數(shù)據(jù)見附錄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 續(xù) #例

2、Box.test(tdat1.5,type=Ljung-Box,lag=12) tdat1.2=Data1.2,2a1.2=acf(tdat1.2)模型擬合序列2.5 續(xù)選擇合適的ARMA#例 acf(tdat1.5)pacf(tdat1.5)根據(jù)自相關(guān)系數(shù)圖和偏自相關(guān)系數(shù)圖可以判斷為12 月平均每頭奶#年例 #2.2繪制1962 年1 月至1975）模型牛產(chǎn)奶量序列時序圖（數(shù)據(jù)見附錄1.3）（1ARData1.3=read.csv(C:UsersAdministratorDesktop口徑的求法在文檔上P81 # 例 2.5續(xù)#P831.3.csv,header=F)附錄然大似陣tdat1.

3、3=as.vector(t(as.matrix(Data1.3)1:168#矩轉(zhuǎn)arima(tdat1.5,order=c(1,0,0),method=ML)#極 置轉(zhuǎn)向量估計ar1=arima(tdat1.5,order=c(1,0,0),method=ML) plot(tdat1.3,type=l)summary(ar1)例 2.2 續(xù) #ev=ar1$residualsacf(tdat1.3)#把字去掉acf(ev) pacf(tdat1.3)pacf(ev) 年北京市每年最高氣溫序19982.3# 例繪制 1949 列時序圖Data1.4=read.csv(C:UsersAdminis

4、tratorDesktop #參數(shù)的顯著性檢驗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# 例續(xù)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 續(xù) Box.test(tdat1.4,type=Ljung-Box,lag=6)P94 續(xù)預(yù)測及置信區(qū)間 #例 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) 個服從標準正態(tài)分布的白噪聲#隨機產(chǎn)生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)#例續(xù) 2.4lines(U,col=lue,lty=dashed)a2.4=acf(Data2.4)lines(L,col=lue,lty=dashed)2.4 例 #續(xù) Box.test(Data2.4,type=Ljung-Box,lag=6)Box.test(Data2.4,type=Ljung

7、-Box,lag=12)續(xù)例例3.53.53.1.1 # 例plot.ts(arima.sim(n=100,list(ar=0.8) #方法一#年北京市城鄉(xiāng)居民定期儲蓄1998 1950 對 2.5例方法二#x0=runif(1)所占比例序列的平穩(wěn)性與純隨機性進行檢驗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#例求平穩(wěn)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)自相關(guān)系數(shù)圖截尾和偏自相關(guān)系模型例題 #3.6 MAfor(i in 3:length(x)數(shù)圖拖尾 #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 必需為數(shù)值for(i in 1:1000)年 1985 模型擬合 1880 #例 3.9 選擇合適的ARMAxi=rnormi-4/5*rnormi -1+16/25*rnormi -2全球氣溫改變差值差分序列plot(x,type=l)沒有數(shù)據(jù) #acf(x)矩估計例 3.12#例 3.10例 3.11pacf(x)次化學反應(yīng)的過程數(shù)xi = rnormi- 4/5 *

14、rnormi - 1 + 16/25 * #錯誤于 703.13 對等時間間隔的連續(xù) #例 rnormi - 2 : 據(jù)進行擬合Data1.8=read.csv(C:UsersAdministratorDesktop #更換參數(shù)長度為零 1.8.csv,header=F)3.6#例續(xù) 64 根據(jù)書上頁來判斷附錄tdat1.8=as.vector(t(as.matrix(Data1.8)1:70plot(tdat1.8,type=o)，（合 3.7#例擬 ARMA11型）模并直觀觀察該模型自相關(guān)，x(t) -0.5x(t -1)=u(t) -0.8*(u -1)）模型 3）例 3.16AR（ 系

15、數(shù)和偏自相關(guān)系數(shù)的拖尾性。 3.15AR例 3.14AR（2）例（ 3# #的預(yù)測法一：x0=runif(1)。 #如果考得話就先。 。 x=rep(0,1000)x1=0.5*x0+rnorm(1)-0.8*rnorm(1)4.1 線性擬合消費支出數(shù)據(jù)#例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) # 法二 #返回擬合參數(shù)的統(tǒng)計量 x=arima.sim(n=1000,list(ar=0.5,ma= -0.8)返回被估計的系數(shù)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)轉(zhuǎn)陣tdat1.9=as.vector(t(as.matrix(Data1.9)1:130#矩擬合圖和原來的圖畫在一起lines(that,type=l)置轉(zhuǎn)向量#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(xiàn)警告信息： :用 # 此外 x=Data1.2#In Ops.factor(left, right) : *對因子沒有意義plot(x,type=o)dx=diff(x,2)指數(shù)平滑法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#例季節(jié)效應(yīng)分析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 #求季節(jié)因子x=ts(Data1.12,2)seaindplot(seaind,type=b) #一階差分季節(jié)指數(shù)圖

22、#dx=diff(x,lag=1,differences=1)plot(Data1.12 -1,1,dx,type=o) 消除季節(jié)因子的影響 noseandat=tdat1.11/seaind # 二階差分 消除季節(jié)因子之 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)平均每頭奶牛產(chǎn)奶量#例 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過差分#沒做好做預(yù)測# 例 5.6續(xù)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)個數(shù)據(jù)隨機抽取 100lines(x,1,x,2-1.96*sqrt(61.95)xt=cumsum(r) # 由隨機游走公式的出的模型公式?jīng)]畫出來 14#圖 5plot(xt,type=l)#隨機游走的圖形5.6

26、續(xù) #dx=diff(xt) #做一階差分例 m3=arima(x,2,order=c(0,1,1),method=ML)plot(dx,type=l) #一階差分后的圖形#m=mean(dx)均值p-171續(xù) 方差sd=var(dx) # #例5.6plot(x,xlim=c(1950,1990),ylim=c(0,300),type=o)用統(tǒng)計量檢驗隨機#Box.test(dx,lag=12,type=Ljung)的函數(shù),data=x) # 變量為時間tm1=lm(農(nóng)業(yè)性年份？模型口徑不會算summary(m1)# 自相關(guān)圖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 例 #對中國農(nóng)業(yè)實際國民收入指數(shù)進行建模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 統(tǒng)計量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=序列具有很強的相關(guān)信息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 疏系數(shù)模型婦女參數(shù)估計與檢驗#例 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 乘積季節(jié)模型#例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,