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

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

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

3、", xlab="x 軸 ",ylab='y 軸 ')有關(guān)作圖命令plot 的詳細(xì)介紹可以在R 中輸入 help(plot)數(shù)學(xué)函數(shù)abs，sqrt ：絕對值， 平方根log, log10, log2 , exp ：對數(shù)與指數(shù)函數(shù)sin，cos，tan，asin，acos，atan， atan2：三角函數(shù)sinh， cosh， tanh ， asinh， acosh， atanh：雙曲函數(shù)簡單統(tǒng)計量sum, mean, var, sd, min, max, range, median, IQR（四分位間距） 等為統(tǒng)計量， sort ，order ，

4、rank與排序有關(guān)，其它還有 ave， fivenum ， mad， quantile ， stem 等。下面我們看一看S的統(tǒng)計功能 :> marks <- c(10, 6, 4, 7, 8)> mean(marks)> sd(marks)> min(marks)> max(marks)第一個語句輸入若干數(shù)據(jù)到一個向量，函 c()用來把數(shù)據(jù)組合為一個向量。后面用了幾個函數(shù)來計算數(shù)據(jù)的均值、標(biāo)準(zhǔn)差、最小值、最大值?？梢园讶舾尚忻畋４嬖谝粋€文本文件中，然后用source 函數(shù)來運(yùn)行整個文件：> source("C:/l.R")注意字符

5、串中的反斜杠。例 :計算 6, 4, 7, 8， 10 的均值和標(biāo)準(zhǔn)差，把若干行命令保存在一個文本文件（比如C:1.R）中，然后用source 函數(shù)來運(yùn)行整個文件。a<- c(10, 6, 4, 7, 8)b<-mean(a)c<-sd(a)source("C:/1.R")時間序列數(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, frequen

6、cy = 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ù)輸出writewrite.tablewrite.csv第二講1. 繪制時序圖、自相關(guān)圖例題 2.1d=scan("sha.csv")sha=ts(d,start=1964,freq=1)plot.ts(sha)#繪制時序圖acf(sha,22)#繪制自相關(guān)圖，滯后期數(shù)22pacf(sha,22)#繪制偏自相關(guān)圖，滯后期數(shù)22corr=acf(sh

7、a,22)#保存相關(guān)系數(shù)cov=acf(sha,22,type = "covariance")#保存協(xié)方差圖的保存，單擊選中圖，在菜單欄選中“文件”，再選“另存為” 。同時顯示多個圖：用x11()命令生成一個空白圖，再輸入作圖命令。2. 同時繪制兩組數(shù)據(jù)的時序圖d=read.csv("double.csv",header=F)double=ts(d,start=1964,freq=1)plot(double, plot.type = "multiple")#兩組數(shù)據(jù)兩個圖plot(double, plot.type = "s

8、ingle")#兩組數(shù)據(jù)一個圖plot(double, plot.type = "single",col=c("red","green"),lty=c(1,2) # 設(shè)置每組數(shù)據(jù)圖的顏色、曲線類型 )3.產(chǎn)生服從正態(tài)分布的隨機(jī)觀察值例題 2.4 隨機(jī)產(chǎn)生1000 白噪聲序列觀察值d=rnorm(1000,0,1)#個數(shù) 1000 均值 0 方差 1plot.ts(d)4.純隨機(jī)性檢驗例題 2.3 續(xù)d=scan("temp.csv")temp=ts(d,freq=1,start=c(1949)Box.te

9、st(temp, type="Ljung-Box",lag=6)5.差分計算x=1:10y=diff(x)k 步差分k xtxtxt k 加入?yún)?shù) lag=k如計算 x 的 3 步差分為y=diff(x, lag = 3)p 階差分p xtp 1xtp 1 xt 1 加入?yún)?shù) differences = p如2階差分2 xtxtxt 1y=diff(x,differences = 2)第三講例題 3.1plot.ts(arima.sim(n = 100, list(ar = 0.8)#模擬 AR(1)模型，并作時序圖。plot.ts(arima.sim(n = 100, l

10、ist(ar = -1.1)#非平穩(wěn)，無法得到時序圖。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.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

11、(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.5#導(dǎo)入數(shù)據(jù)prop=ts(d,start=1950,freq=1)#轉(zhuǎn)化為時間序列數(shù)據(jù)plot(prop)#作時序圖acf(prop,12)#作自相關(guān)圖，拖尾pacf(prop,12)#作偏自相關(guān)圖，1 階截尾Box.test(prop, type="Ljung-Box",lag=

12、6)#純隨機(jī)性檢驗,p 值小于 5%，序列為非白噪聲Box.test(prop, type="Ljung-Box",lag=12)arima(prop, order = c(1,0,0),method="ML")#用 AR(1)模型擬合，如參數(shù) method="CSS"，估計方法為條件最小二乘法，用條件最小二乘法時，不顯示 AIC。arima(prop, order = c(1,0,0),method="ML", include.mean = F) # 用 AR(1)模型擬合， 不含截距項。 tsdiag(arim

13、a(prop, order = c(1,0,0),method="ML")#對估計進(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)# 對殘差進(jìn)行純隨機(jī)性檢驗predict(arima(prop, order = c(1,0,0),

14、n.ahead =5) #預(yù)測未來 5 期 prop.fore = predict(arima(prop, order = c(1,0,0), n.ahead =5) #將未來 5 期預(yù)測值保存在 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)#作時序圖，含預(yù)測。lines(U, col="blue", lty="dashed"

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

16、, 0, 0), method = "ML")Coefficients:ar1intercept0.691481.5509s.e.0.09891.7453sigma2 estimated as 15.51:log likelihood = -137.02,aic = 280.05注意： intercept 下面的 81.5509 是均值，而不是截距！雖然intercept 是截距的意思，這里如果用 mean 會更好。（ the mean and the intercept are the same only when there is no AR term，均值和截距是相同

17、的，只有在沒有AR項的時候）如果想得到截距，利用公式計算。int=(1-0.6914)*81.5509= 25.16661 。課本 P81 的例 2.5 續(xù)中的截距 25.17 是正確的。第二點：如需計算參數(shù)的 t 統(tǒng)計量值和 p 值，利用下面的公式。ar 的 t 統(tǒng)計量值 =0.6914/0.0989= 6.9909(注：數(shù)值與課本略有不同，因為課本用sas 算的 se= 0.1029， R 計算的 se=0.0989)p 值 =pt(6.9909,df=48,lower.tail = F)*2pt() 為求 t 分布求 p 值的函數(shù)， 6.99為 t 統(tǒng)計量的絕對值，df 為自由度 =數(shù)據(jù)

18、個數(shù) -參數(shù)個數(shù)，lower.tail = F 表示所求 p 值為 PT > t，如不加入這個參數(shù)表示所求p 值為 PT <=t。乘 2 表示 p 值是雙側(cè)的（課本上的p 值由 sas 算出，是雙側(cè)的）均值的 t 統(tǒng)計量值和 p 值同理。在時間序列中對參數(shù)顯著性的要求與回歸模型不同，我們更多的是考察模型整體的好壞，而不是參數(shù)。所以， R中的 arima 擬合結(jié)果中沒有給出參數(shù)的t 統(tǒng)計量值和 p 值，如果題目沒有特別要求，一般不需要手動計算。第三點：修正第三講中的錯誤：例 2.5 中，我們用下面的語句對擬合arima 模型之后的殘差進(jìn)行了LB 檢驗：a=arima(prop, or

19、der = c(1,0,0),method="ML")r=a$residualsa=arima(prop, order = c(1,0,0),method="ML")r=a$residuals#用r 來保存殘差Box.test(r,type="Ljung-Box",lag=6)#對殘差進(jìn)行純隨機(jī)性檢驗最后一句不完整，需要加上參數(shù)fitdf=1 ，修改為Box.test(r,type="Ljung-Box",lag=6 ， fitdf=1)fitdf 表示 p+q， number of degrees of free

20、dom to be subtracted if x is a series of residuals驗的序列是殘差到時候，需要加上命令fitdf ，表示減去的自由度。運(yùn)行 Box.test(r,type="Ljung-Box",lag=6,fitdf=1) 后，顯示的結(jié)果：Box.test(r,type="Ljung-Box",lag=6,fitdf=1)Box-Ljung test，當(dāng)檢data:rX-squared = 5.8661, df = 5, p-value = 0.3195“ df = 5表示”自由度為5，由于參數(shù)lag=6，所以是滯后6

21、期的檢驗。第四講# 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)ts.plot(a)t=1:12tlm1=lm(at)summary(lm1)# 返回擬合參數(shù)的統(tǒng)計量coef(lm1)#返回被估計的系數(shù)fitted(lm1)#返回模擬值residuals(lm1) # 返回殘差值fit1=as.ts(fitted(lm1)ts.plot(a);lines(fit1,col="red") #擬合圖

22、#eg1cs=ts(scan("eg1.txt",sep=",")csts.plot(cs)t=1:40lm2=lm(cst)summary(lm2)# 返回擬合參數(shù)的統(tǒng)計量coef(lm2)#返回被估計的系數(shù)fit2=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.8

23、1,528.23,1040.27,2064.25,4113.73,8212.21,16405.95)x2plot(t,x2)m1=nls(x2a*t+bt,start=list(a=0.1,b=1.1),trace=T)summary(m1)# 返回擬合參數(shù)的統(tǒng)計量coef(m1)#返回被估計的系數(shù)fitted(m1)#返回模擬值residuals(m1) # 返回殘差值plot(t,x2);lines(t,fitted(m1) #擬合圖#讀取 excel 中讀取文件，逗號分隔符a=read.csv("example4_2.csv",header=TRUE)t=a$tx=a

24、$xxts.plot(x)m2=nls(xa*t+bt,start=list(a=0.1,b=1.1),trace=T)summary(m2)# 返回擬合參數(shù)的統(tǒng)計量coef(m2)#返回被估計的系數(shù)fitted(m2)#返回模擬值residuals(m2) # 返回殘差值plot(t,x);lines(t,fitted(m2) #擬合圖#eg2I<-scan("eg2.txt")Ix=ts(data=I,start=c(1991,1),f=12) # 化為時間序列xplot.ts(x)t=1:130t2=t2m3=lm(xt+t2)coef(m3)#返回被估計的系數(shù)

25、summary(m3)# 返回擬合參數(shù)的統(tǒng)計量#去不顯著的自變量，再次模擬m4=lm(xt2)coef(m4)summary(m4)m2=fitted(m4)#返回被估計的系數(shù)# 返回擬合參數(shù)的統(tǒng)計量#返回模擬值y=ts(data=m2,start=c(1991,1),f=12)yts.plot(x);lines(y)#平滑法#簡單移動平均法x=c(5,5.4,5.8,6.2)xy=filter(x,rep(1 /4,4),sides=1)y#指數(shù)平滑for(i in 1:3)x1=x1xi+1=0.25*xi+1+0.75*xi#HoltWinters Filtera=ts(read.csv

26、("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(m)plot(fitted(m)#綜合cs=ts(read.csv("eg3.csv",header=F),start=c(1993,1),f=12) #讀取數(shù)據(jù) csts.plot(cs) # 繪制時序圖cs.sea1=rep(0,12)cs.sea1for(i in 1:12)for(j

27、in 1:8)cs.sea1i=cs.sea1i+csi+12*(j-1)cs.sea=(cs.sea1/8)/(mean(cs)cs.seacs.sea2=rep(cs.sea,8)cs.sea2x=cs/cs.sea2xplot(x)t=1:96m1=lm(xt)coef(m1)summary(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)#白噪聲檢驗第五講#回顧#例 5.1sha=ts(s

28、can("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.3milk=ts(scan("milk.txt"),

29、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(dm1)#例 5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)#擬合 ARIMA 模型a=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=arim

30、a(a,c(1,1,0),method="ML")predict(fit1,5)#incom=ts(read.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.plot(dincom)acf(dincom,lag=18)#自相關(guān)圖Box.test(dincom,type="Ljung-Box",lag=6) # 白噪聲檢驗Box.test(dincom,type="Ljung-Box",lag=12

31、)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),xreg=1:length(incom),method="CSS")#見Box.test(fit2$resid,lag=6,type="Ljung-Box",fitdf=1) fore=predict(fit2,10,newxreg=(leng

32、th(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,type="

33、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")summary(fit4)#例 5.9ue=ts(scan("unemployment.txt"),start=1962,f=4) #讀取數(shù)據(jù)par(mfrow=c(2,2) # 繪制時序圖ts.plot(ue)#差分due=

34、diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)#平穩(wě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=F,method="ML")

35、#參數(shù)估計與檢驗（加載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)#殘差白噪聲檢驗Box.test(arma$resid,6,fitdf=2,type="Ljung")#擬合ft=ts(

36、fitted(fit3),start=1963.25,f=4)dft=ts(rep(0,115),start=1963.25,f=4)for(i in 1:115)dfti=fti+duei+uei+4ts.plot(ue);lines(dft,col="red")#例 5.10 乘積季節(jié)模型wue=ts(scan("wue.txt"),start=1948,f=12)arima(wue,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12),include.mean=F,method="C

37、SS")#擬合 Auto-Regressive 模型eg1=ts(scan("582.txt")ts.plot(eg1)#因變量關(guān)于時間的回歸模型fit.gls=gls(eg1-1+time(eg1),correlation=corARMA(p=1),method="ML")#see the nlme packagesummary(fit.gls2)#the results#延遲因變量回歸模型leg1=lag(eg1,-1)y=cbind(eg1,leg1)fit=arima(y,1,c(0,0,0),xreg=y,2,include.mean

38、=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.3milk=ts(scan

39、("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(dm1)#例 5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)#擬合 ARIMA 模型#上機(jī)指導(dǎo)a=ts(scan("581.txt")par(mfrow=c(2,2)ts.plot(a)da=diff(a)ts.plot(da)acf(da,20)pacf(d

40、a,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)# 截距項不顯著，故舍去fit3=arima(a,c(1,1,0),method="ML")predict(fit3,5)#例 5.8incom=ts(rea

41、d.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.plot(dincom)acf(dincom,lag=18)#自相關(guān)圖Box.test(dincom,type="Ljung-Box",lag=6) # 白噪聲檢驗pacf(dincom,lag=18)fit=arima(incom,order=c(0,1,1),xreg=1:length(incom),method="CSS")#見/Rissues.htmAuto

42、corTest(fit$resid)#加載FinTS 包fore=predict(fit,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,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),

43、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),fixed=c(NA,0,0,NA),include.mean=F,data=dw,method="CSS")# 加 載fArma包 ，檢驗參數(shù)summary(fit4)#例 5.9#讀取數(shù)據(jù)ue=ts(scan("unemp

44、loyment.txt"),start=1962,f=4)#繪制時序圖par(mfrow=c(2,2)ts.plot(ue)#差分due=diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)#平穩(wě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=F,method="ML")#參數(shù)估計與檢驗（加載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,meth