R語言常用上機(jī)命令分功能整理——時(shí)間序列分析為主

1、R語言常用上機(jī)命令分功能整理一一時(shí)間序列分析為主第一講應(yīng)用實(shí)例?R的基本界面是一個(gè)交互式命令窗口，命令提示符是一個(gè)大于號(hào)，命令的結(jié)果馬上顯示在命令下面。?S命令主要有兩種形式：表達(dá)式或賦值運(yùn)算(用W或者表示)。在命令提示符后鍵入一個(gè)表達(dá)式表示計(jì)算此表達(dá)式并顯示結(jié)果。賦值運(yùn)算把賦值號(hào)右邊的值計(jì)算出來賦給左邊的變量。?可以用向上光標(biāo)鍵來找回以前運(yùn)行的命令再次運(yùn)行或修改后再運(yùn)行。?S是區(qū)分大小寫的，所以x和X是不同的名字。我們用一些例子來看R軟件的特點(diǎn)。假設(shè)我們已經(jīng)進(jìn)入了R的交互式窗口。如果沒有打開的圖形窗口，在R中，用：x11()可以打開一個(gè)作圖窗口。然后，輸入以下語句：x1=0:100x2=x1

2、*2*pi/100y=sin(x2)plot(x2,y,type="l")這些語句可以繪制正弦曲線圖。其中，“建賦值運(yùn)算符。0:100表示一個(gè)從0到100的等差數(shù)列向量。第二個(gè)語句可以看出，我們可以對(duì)向量直接進(jìn)行四則運(yùn)算，計(jì)算得到的x2是向量x1的所有元素乘以常數(shù)2*pi/100的結(jié)果。從第三個(gè)語句可看到函數(shù)可以以向量為輸入，并可以輸出一個(gè)向量，結(jié)果向量y的每一個(gè)分量是自變量x2的每一個(gè)分量的正弦函數(shù)值。plot（x2,y,type="l",main="畫圖練習(xí)",sub="好好練",xlab="x軸&q

3、uot;,ylab='y軸'）有關(guān)作圖命令plot的詳細(xì)介紹可以在R中輸入help(plot)數(shù)學(xué)函數(shù)abs,sqrt:絕對(duì)值，平方根log,log10,log2,exp：對(duì)數(shù)與指數(shù)函數(shù)sin,cos,tan,asin,acos,atan,atan2：三角函數(shù)sinh,cosh,tanh,asinh,acosh,atanh：雙曲函數(shù)簡單統(tǒng)計(jì)量sum,mean,var,sd,min,max,range,median,IQR(四分位間品巨)等為統(tǒng)計(jì)量，sort,order,rank與排序有關(guān)，其它還有ave,fivenum,mad,quantile,stem等。下面我們看一看S的統(tǒng)

4、計(jì)功能：> marks<-c(10,6,4,7,8)> mean(marks)> sd(marks)> min(marks)> max(marks)第一個(gè)語句輸入若干數(shù)據(jù)到一個(gè)向量，函c()用來把數(shù)據(jù)組合為一個(gè)向量。后面用了幾個(gè)函數(shù)來計(jì)算數(shù)據(jù)的均值、標(biāo)準(zhǔn)差、最小值、最大值?？梢园讶舾尚忻畋４嬖谝粋€(gè)文本文件中，然后用source函數(shù)來運(yùn)行整個(gè)文件：> source("C:/l.R")注意字符串中的反斜杠。例:計(jì)算6,4,7,8,10的均值和標(biāo)準(zhǔn)差，把若干行命令保存在一個(gè)文本文件(比如C:1.R)中，然后用source函數(shù)來運(yùn)行整個(gè)文

5、件。a<-c(10,6,4,7,8)b<-mean(a)c<-sd(a)source("C:/1.R")時(shí)間序列數(shù)據(jù)的輸入使用函數(shù)tsts(1:10,frequency=4,start=c(1959,2)print(ts(1:10,frequency=7,start=c(12,2),calendar=TRUE)a<-ts(1:10,frequency=4,start=c(1959,2)plot(a)將外部數(shù)據(jù)讀入Rread.csv默認(rèn)header=TRUE,也就是第一行是標(biāo)簽，不是數(shù)據(jù)。read.table默認(rèn)header=FALSE將R中的數(shù)據(jù)輸出

6、writewrite.tablewrite.csv第二講1 .繪制時(shí)序圖、自相關(guān)圖例題2.1d=scan("sha.csv")sha=ts(d,start=1964,freq=1)plot.ts(sha)#繪制時(shí)序圖acf(sha,22)#繪制自相關(guān)圖，滯后期數(shù)22pacf(sha,22)#繪制偏自相關(guān)圖，滯后期數(shù)22corr=acf(sha,22)#保存相關(guān)系數(shù)cov=acf(sha,22,type="covariance")琳存協(xié)方差圖的保存，單擊選中圖，在菜單欄選中“文件”，再選“另存為”。同時(shí)顯示多個(gè)圖：用x11()命令生成一個(gè)空白圖，再輸入作圖

7、命令。2 .同時(shí)繪制兩組數(shù)據(jù)的時(shí)序圖d=read.csv("double.csv",header=F)double=ts(d,start=1964,freq=1)plot(double,plot.type="multiple")#兩組數(shù)據(jù)兩個(gè)圖plot(double,plot.type="single")#兩組數(shù)據(jù)一個(gè)圖plot(double,plot.type="single",col=c("red","green"),lty=c(1,2)#設(shè)置每組數(shù)據(jù)圖的顏色、曲線類型)

8、3 .產(chǎn)生服從正態(tài)分布的隨機(jī)觀察值例題2.4隨機(jī)產(chǎn)生1000白噪聲序列觀察值d=rnorm(1000,0,1)#個(gè)數(shù)1000均值0方差1plot.ts(d)4 .純隨機(jī)性檢驗(yàn)例題2.3續(xù)d=scan("temp.csv")temp=ts(d,freq=1,start=c(1949)Box.test(temp,type="Ljung-Box",lag=6)5 .差分計(jì)算x=1:10y=diff(x)k步差分k"xtxtk加入?yún)?shù)lag=k如計(jì)算x的3步差分為y=diff(x,lag=3)pp1p1p階差分xtxtxt1加入?yún)?shù)difference

9、s=p如2階差分僅XtXt1y=diff(x,differences=2)第三講例題3.1plot.ts(arima.sim(n=100,list(ar=0.8)# 模才AAR(1)模型，并作時(shí)序圖。plot.ts(arima.sim(n=100,list(ar=-1.1)# 非平穩(wěn)，無法得到時(shí)序圖。plot.ts(arima.sim(n=100,list(ar=c(1,-0.5)plot.ts(arima.sim(n=100,list(ar=c(1,0.5)例題3.5acf(arima.sim(n=100,list(ar=0.8)acf(arima.sim(n=100,list(ar=-1.

10、1)acf(arima.sim(n=100,list(ar=c(1,-0.5)acf(arima.sim(n=100,list(ar=c(1,0.5)例題3.7arima.sim(n=1000,list(ar=0.5,ma=-0.8)acf(arima.sim(n=1000,list(ar=0.5,ma=-0.8),20)pacf(arima.sim(n=1000,list(ar=0.5,ma=-0.8),20)例題2.5d=scan("a1.5.txt")#導(dǎo)入數(shù)據(jù)prop=ts(d,start=1950,freq=1)#轉(zhuǎn)化為時(shí)間序列數(shù)據(jù)plot(prop)#作時(shí)序圖a

11、cf(prop,12)#作自相關(guān)圖，拖尾pacf(prop,12)#作偏自相關(guān)圖，1階截尾Box.test(prop,type="Ljung-Box",lag=6)# 純隨機(jī)性檢驗(yàn),p值小于5%,序列為非白噪聲Box.test(prop,type="Ljung-Box",lag=12)arima(prop,order=c(1,0,0),method="ML")# 用AR(1膜型擬合，如參數(shù)method="CSS”,估計(jì)方法為條件最小二乘法，用條件最小二乘法時(shí)，不顯示AIC。arima(prop,order=c(1,0,0),

12、method="ML",include.mean=F)#用AR(1)模型擬合，不含截距項(xiàng)。tsdiag(arima(prop,order=c(1,0,0),method="ML")# 對(duì)估計(jì)進(jìn)行診斷，判斷殘差是否為白噪聲summary(arima(prop,order=c(1,0,0),method="ML")a=arima(prop,order=c(1,0,0),method="ML")r=a$residuals#用r來保存殘差Box.test(r,type="Ljung-Box",lag=6

13、)#對(duì)殘差進(jìn)行純隨機(jī)性檢驗(yàn)predict(arima(prop,order=c(1,0,0),n.ahead=5)#預(yù)測(cè)未來5期prop.fore=predict(arima(prop,order=c(1,0,0),n.ahead=5)# 將未來5期預(yù)測(cè)值保存在prop.fore變量中U=prop.fore$pred+1.96*prop.fore$seL=prop.fore$pred-1.96*prop.fore$se#算出95%置信區(qū)間ts.plot(prop,prop.fore$pred,col=1:2)#作時(shí)序圖，含預(yù)測(cè)。lines(U,col="blue",lty=

14、"dashed")lines(L,col="blue",lty="dashed")#在時(shí)序圖中作出95%置信區(qū)間例題3.9d=scan("a1.22.txt")x=diff(d)arima(x,order=c(1,0,1),method="CSS")tsdiag(arima(x,order=c(1,0,1),method="CSS")第一點(diǎn)：對(duì)于第三講中的例2.5,運(yùn)行命令arima(prop,order=c(1,0,0),method="ML")之后，顯

15、示：Call:arima(x=prop,order=c(1,0,0),method="ML")Coefficients:ar1intercept0.691481.5509s.e.0.09891.7453sigmaA2estimatedas15.51:loglikelihood=-137.02,aic=280.05注意：intercept下面的81.5509是均值，而不是截距！雖然intercept是截距的意思，這里如果用mean會(huì)更好。(themeanandtheinterceptarethesameonlywhenthereisnoARterm,均值和截距是相同的，只有在

16、沒有AR項(xiàng)的時(shí)候)如果想得到截距，利用公式計(jì)算。int=(1-0.6914)*81.5509=25.16661。課本P81的例2.5續(xù)中的截距25.17是正確的。第二點(diǎn)：如需計(jì)算參數(shù)的t統(tǒng)計(jì)量值和p值，利用下面的公式。ar的t統(tǒng)計(jì)量值=0.6914/0.0989=6.9909(注：數(shù)值與課本略有不同，因?yàn)檎n本用sas算的se=0.1029,R計(jì)算的se=0.0989)p值=pt(6.9909,df=48,lower.tail=F)*2pt()為求t分布求p值的函數(shù)，6.99為t統(tǒng)計(jì)量的絕對(duì)值，df為自由度=數(shù)據(jù)個(gè)數(shù)-參數(shù)個(gè)數(shù)，lower.tail=F表示所求p值為PT>t,如不加入這個(gè)參

17、數(shù)表示所求p值為PT<=t。乘2表示p值是雙側(cè)的(課本上的p值由sas算出，是雙側(cè)的)均值的t統(tǒng)計(jì)量值和p值同理。在時(shí)間序列中對(duì)參數(shù)顯著性的要求與回歸模型不同，我們更多的是考察模型整體的好壞，而不是參數(shù)。所以，R中的arima擬合結(jié)果中沒有給出參數(shù)的t統(tǒng)計(jì)量值和p值，如果題目沒有特別要求，一般不需要手動(dòng)計(jì)算。第三點(diǎn)：修正第三講中的錯(cuò)誤：例2.5中，我們用下面的語句對(duì)擬合arima模型之后的殘差進(jìn)行了LB檢驗(yàn)：a=arima(prop,order=c(1,0,0),method="ML")r=a$residualsa=arima(prop,order=c(1,0,0),

18、method="ML")r=a$residuals#用r來保存殘差Box.test(r,type="Ljung-Box",lag=6)#對(duì)殘差進(jìn)行純隨機(jī)性檢驗(yàn)最后一句不完整，需要加上參數(shù)fitdf=1,修改為Box.test(r,type="Ljung-Box",lag=6,fitdf=1)fitdf表示p+q,numberofdegreesoffreedomtobesubtractedifxisaseriesofresiduals,當(dāng)檢驗(yàn)的序列是殘差到時(shí)候，需要加上命令fitdf,表示減去的自由度。運(yùn)行Box.test(r,type

19、="Ljung-Box",lag=6,fitdf=1)后，顯示的結(jié)果：Box.test(r,type="Ljung-Box",lag=6,fitdf=1)Box-Ljungtestdata:rX-squared=5.8661,df=5,p-value=0.3195“df=裱市自由度為5,由于參數(shù)lag=6,所以是滯后6期的檢驗(yàn)。第四講#example4_1擬合線性模型x1=c(12.79,14.02,12.92,18.27,21.22,18.81,25.73,26.27,26.75,28.73,31.71,33.95)a=as.ts(x1)is.ts(a

20、)ts.plot(a)t=1:12tlm1=lm(at)summary(lm1)#返回?cái)M合參數(shù)的統(tǒng)計(jì)量coef(lm1)#返回被估計(jì)的系數(shù)fitted(lm1)#返回模擬值residuals(lm1)#返回殘差值fit1=as.ts(fitted(lm1)ts.plot(a);lines(fit1,col="red")#擬合圖#eg1cs=ts(scan("eg1.txt",sep=",")csts.plot(cs)t=1:40lm2=lm(cst)summary(lm2)#返回?cái)M合參數(shù)的統(tǒng)計(jì)量coef(lm2)#返回被估計(jì)的系數(shù)fi

21、t2=as.ts(fitted(lm2)#返回模擬值residuals(lm2)#返回殘差值ts.plot(cs);lines(fit2,col="red")#擬合圖#example4_2擬合非線性模型t=1:14x2=c(1.85,7.48,14.29,23.02,37.42,74.27,140.72,265.81,528.23,1040.27,2064.25,4113.73,8212.21,16405.95)x2plot(t,x2)m1=nls(x2-a*t+bAt,start=list(a=0.1,b=1.1),trace=T)summary(ml)#返回?cái)M合參數(shù)的統(tǒng)

22、計(jì)量coef(m1)#返回被估計(jì)的系數(shù)fitted(ml)#返回模擬值residuals(ml)#返回殘差值plot(t,x2);lines(t,fitted(m1)#擬合圖#讀取excel中讀取文件，逗號(hào)分隔符a=read.csv("example4_2.csv",header=TRUE)t=a$tx=a$xxts.plot(x)m2=nls(xa*t+bAt,start=list(a=0.1,b=1.1),trace=T)summary(m2)#返回?cái)M合參數(shù)的統(tǒng)計(jì)量coef(m2)#返回被估計(jì)的系數(shù)fitted(m2)#返回模擬值residuals(m2)#返回殘差值p

23、lot(t,x);lines(t,fitted(m2)#擬合圖#eg2I<-scan("eg2.txt")Ix=ts(data=I,start=c(1991,1),f=12)#化為時(shí)間序列xplot.ts(x)t=1:130t2=tA2m3=lm(xt+t2)coef(m3)#返回被估計(jì)的系數(shù)summary(m3)#返回?cái)M合參數(shù)的統(tǒng)計(jì)量# 去不顯著的自變量，再次模擬m4=lm(xt2)coef(m4)#返回被估計(jì)的系數(shù)summary(m4)#返回?cái)M合參數(shù)的統(tǒng)計(jì)量m2=fitted(m4)#返回模擬值y=ts(data=m2,start=c(1991,1),f=12)y

24、ts.plot(x);lines(y)# 平滑法# 簡單移動(dòng)平均法x=c(5,5.4,5.8,6.2)xy=filter(x,rep(1/4,4),sides=1)y# 指數(shù)平滑for(iin1:3)x1=x1xi+1=0.25*xi+1+0.75*xi#HoltWintersFiltera=ts(read.csv("holt.csv",header=F),start=c(1978,1),f=1)am=HoltWinters(a,alpha=0.15,beta=0.1,gamma=FALSE,l.start=51259,b.start=4325)mfitted(m)plot

25、(m)plot(fitted(m)# 綜合cs=ts(read.csv("eg3.csv",header=F),start=c(1993,1),f=12)#讀取數(shù)據(jù)csts.plot(cs)#繪制時(shí)序圖cs.sea1=rep(0,12)cs.sea1for(iin1:12)for(jin1:8)cs.sea1i=cs.sea1i+csi+12*(j-1)cs.sea=(cs.sea/8)/(mean(cs)cs.seacs.sea2=rep(cs.sea,8)cs.sea2x=cs/cs.sea2xplot(x)t=1:96m1=lm(x-t)coef(m1)summary

26、(m1)m=ts(fitted(m1),start=c(1993,1),f=12)ts.plot(x,type="p");lines(m,col="red")r=residuals(m1)Box.test(r)#白噪聲檢驗(yàn)第五講ffffffffffffffffffffffffffffffffffffffffffffffffTT，kk，kk，kk，TT，TTfjf>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>

27、;f>f>f># 回顧# 例5.1sha=ts(scan("sha.csv"),start=1964,freq=1)ts.plot(sha)diff(sha)par(mfrow=c(2,1)ts.plot(diff(sha)acf(diff(sha)# 例5.2car=ts(read.csv("car.csv",header=F),start=1950,freq=1)carpar(mfrow=c(3,1)ts.plot(car)ts.plot(diff(car)ts.plot(diff(car,differences=2)# 例5.3m

28、ilk=ts(scan("milk.txt"),start=c(1962,1),freq=12)milkpar(mfrow=c(3,1)ts.plot(milk)ts.plot(diff(milk)dm1=diff(diff(milk),lag=12)ts.plot(dm1)acf(dml)# 例5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)ffffffffffffffffffffffffffffffffffffffffffffffffffffff，kk，kk，f>f>f>f>f>f>f>f&g

29、t;f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f># 擬合ARIMA模型#5.8.1a=ts(scan("581.txt")par(mfrow=c(2,2)ts.plot(a)da=diff(a)ts.plot(da)acf(da,20)pacf(da,20)Box.test(da,6)fit1=arima(a,c(1,1,0),method="ML")predict(fit1,5)fffffffffffffff

30、fffffffffffffffffffffffffffffffffffffffffff，kk，f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>incom=ts(read.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.

31、plot(dincom)acf(dincom,lag=18)#自相關(guān)圖Box.test(dincom,type="Ljung-Box",lag=6)#白噪聲檢驗(yàn)Box.test(dincom,type="Ljung-Box",lag=12)Box.test(dincom,type="Ljung-Box",lag=18)pacf(dincom,lag=18)fit1=arima(dincom,order=c(0,0,1),method="CSS")fit2=arima(incom,order=c(0,1,1),xre

32、g=1:length(incom),method="CSS")#見/stoffer/tsa2/Rissues.htmBox.test(fit2$resid,lag=6,type="Ljung-Box",fitdf=1)fore=predict(fit2,10,newxreg=(length(incom)+1):(length(incom)+10)# 疏系數(shù)模型# 例5.8w=ts(read.csv("w.csv"),start=1917,freq=1)w=w,1par(mfrow=c(2,

33、2)ts.plot(w)ts.plot(diff(w)acf(diff(w),lag=18)pacf(diff(w),lag=18)dw=diff(w)fit3=arima(dw,order=c(4,0,0),fixed=c(NA,0,0,NA,0),method="CSS")Box.test(fit3$resid,lag=6,type="Ljung-Box",fitdf=2)Box.test(fit3$resid,lag=12,type="Ljung-Box",fitdf=2)fit4=armaFit(arima(4,0,0),fi

34、xed=c(NA,0,0,NA),include.mean=F,data=dw,method="CSS")summary(fit4)# 例5.9ue=ts(scan("unemployment.txt"),start=1962,f=4)#讀取數(shù)據(jù)par(mfrow=c(2,2)#繪制時(shí)序圖ts.plot(ue)# 差分due=diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)# 平穩(wěn)性檢驗(yàn)acf(ddue,lag=30)pacf(ddue,lag=30)arima(ddue,orde

35、r=c(0,0,0),method="ML")arima(ddue,order=c(4,0,0),method="ML")arma=arima(ddue,order=c(4,0,0),transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="ML")# 參數(shù)估計(jì)與檢驗(yàn)(加載fArma程序包)fit2=armaFit(arima(4,0,0),include.mean=F,data=ddue,method="ML")summary(fit2)fit3=arm

36、aFit(arima(4,0,0),data=ddue,transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="CSS")summary(fit3)# 殘差白噪聲檢驗(yàn)Box.test(arma$resid,6,fitdf=2,type="Ljung")# 擬合ft=ts(fitted(fit3),start=1963.25,f=4)dft=ts(rep(0,115),start=1963.25,f=4)for(iin1:115)dfti=fti+duei+uei+4ts.plot(ue);li

37、nes(dft,col="red")ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff“TTk，“TTk，TTk，kk，k，.，ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff# 例5.10乘積季節(jié)模型wue=ts(scan("wue.txt"),start=1948,f=12)arima(wue,order=c(1,1,1),seasonal=list(or

38、der=c(0,1,1),period=12),include.mean=F,method="CSS")ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff，kk，kk，”mrHIIIIIIffffffffffffffffffffffffffffffffffffffffffffffffffffffffffifff# 擬合Auto-Regressive模型eg1=ts(scan("582.txt")ts.plot(eg1)# 因變量關(guān)于時(shí)間的回歸模型fit.gls=gl

39、s(eg1-1+time(eg1),correlation=corARMA(p=1),method="ML")#seethenlmepackagesummary(fit.gls2)#theresults# 延遲因變量回歸模型leg1=lag(eg1,-1)y=cbind(eg1,leg1)fit=arima(y,1,c(0,0,0),xreg=y,2,include.mean=F)第六講# 回顧# 例5.1sha=ts(scan("sha.csv"),start=1964,freq=1)ts.plot(sha)diff(sha)par(mfrow=c(2

40、,1)ts.plot(diff(sha)acf(diff(sha)# 例5.2car=ts(read.csv("car.csv",header=F),start=1950,freq=1)carpar(mfrow=c(3,1)ts.plot(car)ts.plot(diff(car)ts.plot(diff(car,differences=2)# 例5.3milk=ts(scan("milk.txt"),start=c(1962,1),freq=12)milkpar(mfrow=c(3,1)ts.plot(milk)ts.plot(diff(milk)dm

41、1=diff(diff(milk),lag=12)ts.plot(dm1)acf(dm1)# 例5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff# 擬合ARIMA模型# 上機(jī)指導(dǎo)5.8.1a=ts(scan("581.txt")par(mfrow=c(2,2)ts.plot(a)da=diff(a)ts.plot(da)acf

42、(da,20)pacf(da,20)Box.test(da,6)fit1=arima(a,c(1,1,0),method="ML")predict(fit1,5,newxreg=(length(a)+1):(length(a)+5)fit2=armaFit(arima(1,1,0),data=a,xreg=1:length(a),method="ML")summary(fit1)summary(fit2)#截距項(xiàng)不顯著，故舍去fit3=arima(a,c(1,1,0),method="ML")predict(fit3,5)ffffff

43、fffffffffffffffffffffffffffffffffffffffffffffffffffffjf>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f># 例5.8incom=ts(read.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(i

44、ncom)ts.plot(dincom)acf(dincom,lag=18)#自相關(guān)圖Box.test(dincom,type="Ljung-Box",lag=6)#白噪聲檢驗(yàn)pacf(dincom,lag=18)fit=arima(incom,order=c(0,1,1),xreg=1:length(incom),method="CSS")# 見/stoffer/tsa2/Rissues.htmAutocorTest(fit$resid)#力口載FinTS包fore=predict(fit,10,newxr

45、eg=(length(incom)+1):(length(incom)+10)# 疏系數(shù)模型# 例5.8w=ts(read.csv("w.csv"),start=1917,freq=1)w=w,1par(mfrow=c(2,2)ts.plot(w)ts.plot(diff(w)acf(diff(w),lag=18)pacf(diff(w),lag=18)dw=diff(w)fit3=arima(dw,order=c(4,0,0),fixed=c(NA,0,0,NA,0),method="CSS")Box.test(fit3$resid,lag=6,typ

46、e="Ljung-Box",fitdf=2)Box.test(fit3$resid,lag=12,type="Ljung-Box",fitdf=2)fit4=armaFit(arima(4,0,0),fixed=c(NA,0,0,NA),include.mean=F,data=dw,method="CSS")#力口載fArma包，檢驗(yàn)參數(shù)summary(fit4)# 例5.9# 讀取數(shù)據(jù)ue=ts(scan("unemployment.txt"),start=1962,f=4)# 繪制時(shí)序圖par(mfrow=c(

47、2,2)ts.plot(ue)# 差分due=diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)# 平穩(wěn)性檢驗(yàn)acf(ddue,lag=30)pacf(ddue,lag=30)arima(ddue,order=c(0,0,0),method="ML")arima(ddue,order=c(4,0,0),method="ML")arma=arima(ddue,order=c(4,0,0),transform.pars=F,fixed=c(NA,0,0,NA),include.mean=

48、F,method="ML")# 參數(shù)估計(jì)與檢驗(yàn)(加載fArma程序包)fit2=armaFit(arima(4,0,0),include.mean=F,data=ddue,method="ML")summary(fit2)fit3=armaFit(arima(4,0,0),data=ddue,transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="CSS")summary(fit3)# 殘差白噪聲檢驗(yàn)Box.test(arma$resid,6,fitdf=2,type="Ljung")# 擬合ft=ts(fitted(fit3),start=1963.25,f=4)dft=ts(rep(0,115),start=1963.25,f=4)for(iin1:115)dfti=fti+duei+