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

R语言常用上机命令分功能整理——时间序列分析为主第一讲应用实例•R的基本界面是一个交互式命令窗口，命令提示符是一个大于号，命令的结果马上显示在命令下面。•S命令主要有两种形式：表达式或赋值运算（用’－’或者’=’表示）。在命令提示符后键入一个表达式表示计算此表达式并显示结果。赋值运算把赋值号右边的值计算出来赋给左边的变量。•可以用向上光标键来找回以前运行的命令再次运行或修改后再运行。•S是区分大小写的，所以x和X是不同的名字。我们用一些例子来看R软件的特点。假设我们已经进入了R的交互式窗口。如果没有打开的图形窗口，在R中，用：x11()可以打开一个作图窗口。然后，输入以下语句：x1=0:100x2=x1*2*pi/100y=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=好好练,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与排序有关，其它还有ave，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，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,2))plot(a)将外部数据读入Rread.csv默认header=TRUE，也就是第一行是标签，不是数据。read.table默认header=FALSE将R中的数据输出writewrite.tablewrite.csv第二讲1.绘制时序图、自相关图例题2.1d=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(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)#两组数据一个图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-Box,lag=6)5.差分计算x=1:10y=diff(x)k步差分ktttkxxx加入参数lag=k如计算x的3步差分为y=diff(x,lag=3)p阶差分111ppptttxxx加入参数differences=p如2阶差分21tttxxx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.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,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,type=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))#对估计进行诊断，判断残差是否为白噪声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,order=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%置信区间例题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)Coefficients:ar1intercept0.691481.5509s.e.0.09891.7453sigma^2estimatedas15.51:loglikelihood=-137.02,aic=280.05注意：intercept下面的81.5509是均值，而不是截距！虽然intercept是截距的意思，这里如果用mean会更好。（themeanandtheinterceptarethesameonlywhenthereisnoARterm，均值和截距是相同的，只有在没有AR项的时候）如果想得到截距，利用公式计算。int=(1-0.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值为P[Tt]，如不加入这个参数表示所求p值为P[T=t]。乘2表示p值是双侧的（课本上的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)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，numberofdegreesoffreedomtobesubtractedifxisaseriesofresiduals，当检验的序列是残差到时候，需要加上命令fitdf，表示减去的自由度。运行Box.test(r,type=Ljung-Box,lag=6,fitdf=1)后，显示的结果：Box.test(r,type=Ljung-Box,lag=6,fitdf=1)Box