R语言常用上机命令分功能整理——时间序列分析为主

1、R语言常用上机命令分功能整理时间序列分析为主第一讲应用实例 R的基本界面是一个交互式命令窗口，命令提示符是一个大于号，命令的结果马上显示在命令下面。 S命令主要有两种形式：表达式或赋值运算（用<或者=表示）。在命令提示符后键入一个表达式表示计算此表达式并显示结果。赋值运算把赋值号右边的值计算出来赋给左边的变量。 可以用向上光标键来找回以前运行的命令再次运行或修改后再运行。 S是区分大小写的，所以x和X是不同的名字。我们用一些例子来看R软件的特点。假设我们已经进入了R的交互式窗口。如果没有打开的图形窗口，在R中，用：> x11() 可以打开一个作图窗口。然后，输入以下语句： x1 =

2、 0:100 x2 = x1*2*pi/100 y = sin(x2) plot(x2,y,type="l") 这些语句可以绘制正弦曲线图。其中，“=”是赋值运算符。0:100表示一个从0到100 的等差数列向量。第二个语句可以看出，我们可以对向量直接进行四则运算，计算得到的x2 是向量x1的所有元素乘以常数2*pi/100的结果。从第三个语句可看到函数可以以向量为输入，并可以输出一个向量，结果向量y的每一个分量是自变量x2的每一个分量的正弦函数值。plot(x2,y, type="l",main="画图练习",sub="好

3、好练", xlab="x轴",ylab='y轴')有关作图命令plot的详细介绍可以在R中输入help(plot)数学函数abs，sqrt：绝对值，平方根 log, log10, log2 , exp：对数与指数函数 sin，cos，tan，asin，acos，atan，atan2：三角函数 sinh，cosh，tanh，asinh，acosh，atanh：双曲函数 简单统计量sum, mean, var, sd, min, max, range, median, IQR（四分位间距）等为统计量，sort，order，rank与排序有关，其它还有a

4、ve，fivenum，mad，quantile，stem等。下面我们看一看S的统计功能:> marks <- c(10, 6, 4, 7, 8) > mean(marks) > sd(marks) > min(marks)> max(marks) 第一个语句输入若干数据到一个向量，函c()用来把数据组合为一个向量。后面用了几个函数来计算数据的均值、标准差、最小值、最大值。可以把若干行命令保存在一个文本文件中，然后用source函数来运行整个文件：> source("C:/l.R") 注意字符串中的反斜杠。例:计算6, 4, 7, 8

5、，10的均值和标准差，把若干行命令保存在一个文本文件（比如C:1.R）中，然后用source函数来运行整个文件。a<- c(10, 6, 4, 7, 8) b<-mean(a) c<-sd(a) source("C:/1.R")时间序列数据的输入使用函数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

6、, 2)plot(a)将外部数据读入Rread.csv默认header = TRUE，也就是第一行是标签，不是数据。read.table默认header = FALSE将R中的数据输出write write.table write.csv第二讲1. 绘制时序图、自相关图例题2.1 d=scan("sha.csv")sha=ts(d,start=1964,freq=1)plot.ts(sha) #绘制时序图acf(sha,22) #绘制自相关图，滞后期数22pacf(sha,22) #绘制偏自相关图，滞后期数22corr=acf(sha,22) #保存相关系数cov=acf(

7、sha,22,type = "covariance") #保存协方差图的保存，单击选中图，在菜单栏选中“文件”，再选“另存为”。同时显示多个图：用x11()命令生成一个空白图，再输入作图命令。2. 同时绘制两组数据的时序图d=read.csv("double.csv",header=F)double=ts(d,start=1964,freq=1)plot(double, plot.type = "multiple") #两组数据两个图plot(double, plot.type = "single") #两组数据一

8、个图 plot(double, plot.type = "single",col=c("red","green"),lty=c(1,2) #设置每组数据图的颜色、曲线类型)3.产生服从正态分布的随机观察值例题2.4 随机产生1000白噪声序列观察值d=rnorm(1000,0,1) #个数1000 均值0 方差1plot.ts(d)4.纯随机性检验例题2.3续d=scan("temp.csv")temp=ts(d,freq=1,start=c(1949)Box.test(temp, type="Ljung

9、-Box",lag=6)5.差分计算x=1:10y=diff(x)k步差分 加入参数 lag=k如计算x的3步差分为y=diff(x, lag = 3)p阶差分 加入参数differences = p如2阶差分y=diff(x,differences = 2)第三讲例题3.1plot.ts(arima.sim(n = 100, list(ar = 0.8)#模拟AR(1)模型，并作时序图。plot.ts(arima.sim(n = 100, list(ar = -1.1)#非平稳，无法得到时序图。plot.ts(arima.sim(n = 100, list(ar = c(1,-0.

10、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(n = 1000, list(ar = 0.5, ma = -0.8)acf(arima.sim(n = 1000, list(ar = 0.5

11、, ma = -0.8),20)pacf(arima.sim(n = 1000, list(ar = 0.5, ma = -0.8),20)例题2.5d=scan("a1.5.txt") #导入数据prop=ts(d,start=1950,freq=1) #转化为时间序列数据plot(prop) #作时序图acf(prop,12) #作自相关图，拖尾pacf(prop,12) #作偏自相关图，1阶截尾Box.test(prop, type="Ljung-Box",lag=6) #纯随机性检验,p值小于5%，序列为非白噪声Box.test(prop, ty

12、pe="Ljung-Box",lag=12)arima(prop, order = c(1,0,0),method="ML") #用AR(1)模型拟合，如参数method="CSS"，估计方法为条件最小二乘法，用条件最小二乘法时，不显示AIC。arima(prop, order = c(1,0,0),method="ML", include.mean = F) #用AR(1)模型拟合，不含截距项。tsdiag(arima(prop, order = c(1,0,0),method="ML") #

13、对估计进行诊断，判断残差是否为白噪声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)#对残差进行纯随机性检验predict(arima(prop, order = c(1,0,0), n.ahead =5) #预测未来5期prop.fore = predict(arima(prop, orde

14、r = c(1,0,0), n.ahead =5)#将未来5期预测值保存在prop.fore变量中U = prop.fore$pred + 1.96* prop.fore$seL = prop.fore$pred 1.96* prop.fore$se#算出95%置信区间ts.plot(prop, prop.fore$pred,col=1:2)#作时序图，含预测。lines(U, col="blue", lty="dashed")lines(L, col="blue", lty="dashed")#在时序图中作出95

15、%置信区间例题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")第一点：对于第三讲中的例2.5，运行命令arima(prop, order = c(1,0,0),method="ML")之后，显示：Call:arima(x = prop, order = c(1, 0, 0), method = "ML")Coeff

16、icients: ar1 intercept 0.6914 81.5509s.e. 0.0989 1.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，均值和截距是相同的，只有在没有AR项的时候）如果想得到截距，利用公式计算。int=(1-0.

17、6914)*81.5509= 25.16661。课本P81的例2.5续中的截距25.17是正确的。第二点：如需计算参数的t统计量值和p值，利用下面的公式。ar的t统计量值=0.6914/0.0989= 6.9909(注：数值与课本略有不同，因为课本用sas算的se= 0.1029，R计算的se=0.0989)p值=pt(6.9909,df=48,lower.tail = F)*2pt()为求t分布求p值的函数，6.99为t统计量的绝对值，df为自由度=数据个数-参数个数，lower.tail = F表示所求p值为PT > t，如不加入这个参数表示所求p值为PT <=t。乘2表示p值

18、是双侧的（课本上的p值由sas算出，是双侧的）均值的t统计量值和p值同理。在时间序列中对参数显著性的要求与回归模型不同，我们更多的是考察模型整体的好坏，而不是参数。所以，R中的arima拟合结果中没有给出参数的t统计量值和p值，如果题目没有特别要求，一般不需要手动计算。第三点：修正第三讲中的错误：例2.5中，我们用下面的语句对拟合arima模型之后的残差进行了LB检验：a=arima(prop, order = c(1,0,0),method="ML")r=a$residualsa=arima(prop, order = c(1,0,0),method="ML&q

19、uot;)r=a$residuals#用r来保存残差Box.test(r,type="Ljung-Box",lag=6)#对残差进行纯随机性检验最后一句不完整，需要加上参数fitdf=1，修改为Box.test(r,type="Ljung-Box",lag=6，fitdf=1)fitdf表示p+q，number of degrees of freedom to be subtracted if x is a series of residuals，当检验的序列是残差到时候，需要加上命令fitdf，表示减去的自由度。运行Box.test(r,type=&q

20、uot;Ljung-Box",lag=6,fitdf=1)后，显示的结果：Box.test(r,type="Ljung-Box",lag=6,fitdf=1) Box-Ljung testdata: r X-squared = 5.8661, df = 5, p-value = 0.3195“df = 5”表示自由度为5，由于参数lag=6，所以是滞后6期的检验。第四讲# 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=a

21、s.ts(x1)is.ts(a)ts.plot(a)t=1:12tlm1=lm(at)summary(lm1) # 返回拟合参数的统计量coef(lm1) #返回被估计的系数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) # 返回拟合参数

22、的统计量coef(lm2) #返回被估计的系数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.81,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),

23、trace=T)summary(m1) # 返回拟合参数的统计量coef(m1) #返回被估计的系数fitted(m1) #返回模拟值residuals(m1) #返回残差值plot(t,x2);lines(t,fitted(m1) #拟合图#读取excel中读取文件，逗号分隔符 a=read.csv("example4_2.csv",header=TRUE)t=a$tx=a$xxts.plot(x)m2=nls(xa*t+bt,start=list(a=0.1,b=1.1),trace=T)summary(m2) # 返回拟合参数的统计量coef(m2) #返回被估计的系

24、数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) #返回被估计的系数 summary(m3) # 返回拟合参数的统计量#去不显著的自变量 ，再次模拟 m4=lm(xt2) coef(m4) #返回被估计的系数summary(m4) # 返回拟合参数的统计量m2=

25、fitted(m4) #返回模拟值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#指数平滑for(i in 1:3) x1=x1 xi+1=0.25*xi+1+0.75*xi #HoltWinters Filtera=ts(read.csv("holt.csv",header=F),start=c(1978,1),f=1)am=HoltWinters(a,alpha=0.15,beta=

26、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) #读取数据 csts.plot(cs) #绘制时序图 cs.sea1=rep(0,12)cs.sea1for(i in 1:12) for(j in 1:8) cs.sea1i=cs.sea1i+csi+12*(j-1) cs.sea=(cs.sea1/8)/(mean(cs)cs.seacs.sea2=re

27、p(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(scan("sha.csv"),start=1964,freq=1)ts.plot(sha)diff(sha)par(mfrow=c(2,1)

28、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)dm1=dif

29、f(diff(milk),lag=12)ts.plot(dm1)acf(dm1)#例5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)#拟合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)#incom=ts(read.csv("inco

30、m.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.plot(dincom)acf(dincom,lag=18) #自相关图Box.test(dincom,type="Ljung-Box",lag=6) #白噪声检验Box.test(dincom,type="Ljung-Box",lag=12)Box.test(dincom,type="Ljung-Box",lag=18)pacf(dincom,lag=18)fit1=arim

31、a(dincom,order=c(0,0,1),method="CSS")fit2=arima(incom,order=c(0,1,1),xreg=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)

32、#疏系数模型#例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="Ljung-Box",fitdf=2)Box.test(

33、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) #读取数据par(mfrow=c(2,2) #绘制时序图ts.plot(ue)#差分due=diff(ue)ddue=diff(due,lag=4)ts.plo

34、t(ddue)Box.test(ddue,lag=6)#平稳性检验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")#参数估计与检验（加载fArma程序包）fit2=armaFit(a

35、rima(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(fitted(fit3),start=1963.25,f=4)dft=

36、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 乘积季节模型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="CSS")#拟合Auto-Regressive模型eg1=ts

37、(scan("582.txt")ts.plot(eg1)#因变量关于时间的回归模型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=F)第六讲#回顾#例5.1sha=ts(scan("sha.

38、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"),start=c(1962,1),

39、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模型#上机指导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

40、,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.8 incom=ts(read.csv("incom.csv",header=F)

41、,start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.plot(dincom)acf(dincom,lag=18) #自相关图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") #见/stoffer/tsa2/Rissues.htmAuto

42、corTest(fit$resid) #加载FinTS包 fore=predict(fit,10,newxreg=(length(incom)+1):(length(incom)+10)#疏系数模型#例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包 ，检验参数 summary(fit4)#例 5.9#读取数据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)#平稳性检验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")#参数估计与检验（加载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=&