整合R语言常用上机命令分功能整理时间序列分析为主名师精品资料.docx
《整合R语言常用上机命令分功能整理时间序列分析为主名师精品资料.docx》由会员分享，可在线阅读，更多相关《整合R语言常用上机命令分功能整理时间序列分析为主名师精品资料.docx（25页珍藏版）》请在冰豆网上搜索。
整合R语言常用上机命令分功能整理时间序列分析为主名师精品资料
应用实例
•R的基本界面是一个交互式命令窗口，命令提示符是一个大于号，命令的结果马上显示在命令下面。
•S命令主要有两种形式：
表达式或赋值运算（用’<－’或者’=’表示）。
在命令提示符后键入一个表达式表示计算此表达式并显示结果。
赋值运算把赋值号右边的值计算出来赋给左边的变量。
•可以用向上光标键来找回以前运行的命令再次运行或修改后再运行。
•S是区分大小写的，所以x和X是不同的名字。
我们用一些例子来看R软件的特点。
假设我们已经进入了R的交互式窗口。
如果没有打开的图形窗口，在R中，用：
>x11（）可以打开一个作图窗口。
然后，输入以下语句：
x1=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="好好练",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"）
时间序列数据的输入
使用函数ts
ts（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）
将外部数据读入R
read.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）#绘制自相关图，滞后期数22
pacf（sha,22）#绘制偏自相关图，滞后期数22
corr=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方差1
plot.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:
10
y=diff（x）
k步差分
加入参数lag=k
如计算x的3步差分为
y=diff（x,lag=3）
p阶差分
加入参数differences=p
如2阶差分
y=diff（x,differences=2）
第三讲
例题3.1
plot.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.5
acf（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.7
arima.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.5
d=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$se
L=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.9
d=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:
ar1intercept
0.691481.5509
s.e.0.09891.7453
sigma^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）*2
pt（）为求t分布求p值的函数，6.99为t统计量的绝对值，df为自由度=数据个数-参数个数，lower.tail=F表示所求p值为P[T>t]，如不加入这个参数表示所求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$residuals
a=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-Ljungtest
data:
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=as.ts（x1）
is.ts（a）
ts.plot（a）
t=1:
12
t
lm1=lm（a~t）
summary（lm1）#返回拟合参数的统计量
coef（lm1）#返回被估计的系数
fitted（lm1）#返回模拟值
residuals（lm1）#返回残差值
fit1=as.ts（fitted（lm1））
ts.plot（a）;lines（fit1,col="red"）#拟合图
#eg1
cs=ts（scan（"eg1.txt",sep=","））
cs
ts.plot（cs）
t=1:
40
lm2=lm（cs~t）
summary（lm2）#返回拟合参数的统计量
coef（lm2）#返回被估计的系数
fit2=as.ts（fitted（lm2））#返回模拟值
residuals（lm2）#返回残差值
ts.plot（cs）;lines（fit2,col="red"）#拟合图
#example4_2拟合非线性模型
t=1:
14
x2=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）
x2
plot（t,x2）
m1=nls（x2~a*t+b^t,start=list（a=0.1,b=1.1）,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$t
x=a$x
x
ts.plot（x）
m2=nls（x~a*t+b^t,start=list（a=0.1,b=1.1）,trace=T）
summary（m2）#返回拟合参数的统计量
coef（m2）#返回被估计的系数
fitted（m2）#返回模拟值
residuals（m2）#返回残差值
plot（t,x）;lines（t,fitted（m2））#拟合图
#eg2
I<-scan（"eg2.txt"）
I
x=ts（data=I,start=c（1991,1）,f=12）#化为时间序列
x
plot.ts（x）
t=1:
130
t2=t^2
m3=lm（x~t+t2）
coef（m3）#返回被估计的系数
summary（m3）#返回拟合参数的统计量
#去不显著的自变量，再次模拟
m4=lm（x~t2）
coef（m4）#返回被估计的系数
summary（m4）#返回拟合参数的统计量
m2=fitted（m4）#返回模拟值
y=ts（data=m2,start=c（1991,1）,f=12）
y
ts.plot（x）;lines（y）
#平滑法
#简单移动平均法
x=c（5,5.4,5.8,6.2）
x
y=filter（x,rep（1/4,4）,sides=1）
y
#指数平滑
for（iin1:
3）
{
x[1]=x[1]
x[i+1]=0.25*x[i+1]+0.75*x[i]
}
#HoltWintersFilter
a=ts（read.csv（"holt.csv",header=F）,start=c（1978,1）,f=1）
a
m=HoltWinters（a,alpha=0.15,beta=0.1,gamma=FALSE,l.start=51259,b.start=4325）
m
fitted（m）
plot（m）
plot（fitted（m））
#综合
cs=ts（read.csv（"eg3.csv",header=F）,start=c（1993,1）,f=12）#读取数据
cs
ts.plot（cs）#绘制时序图
cs.sea1=rep（0,12）
cs.sea1
for（iin1:
12）{
for（jin1:
8）{
cs.sea1[i]=cs.sea1[i]+cs[i+12*（j-1）]
}
}
cs.sea=（cs.sea1/8）/（mean（cs））
cs.sea
cs.sea2=rep（cs.sea,8）
cs.sea2
x=cs/cs.sea2
x
plot（x）
t=1:
96
m1=lm（x~t）
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.1
sha=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.2
car=ts（read.csv（"car.csv",header=F）,start=1950,freq=1）
car
par（mfrow=c（3,1））
ts.plot（car）
ts.plot（diff（car））
ts.plot（diff（car,differences=2））
#例5.3
milk=ts（scan（"milk.txt"）,start=c（1962,1）,freq=12）
milk
par（mfrow=c（3,1））
ts.plot（milk）
ts.plot（diff（milk））
dm1=diff（diff（milk）,lag=12）
ts.plot（dm1）
acf（dm1）
#例5.5
x=ts（cumsum（rnorm（1000,0,100）））
ts.plot（x）
###########################
#拟合ARIMA模型
#5.8.1
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=arima（a,c（1,1,0）,method="ML"）
predict（fit1,5）
#############################
incom=ts（read.csv（"incom.csv",header=F）,start=1952,freq=1）
incom
ts.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=arima（dincom,order=c（0,0,1）,method="CSS"）
fit2=arima（incom,order=c（0,1,1）,xreg=1:
length（incom）,method="CSS"）#
见http:
//www.stat.pitt.edu/stoffer/tsa2/Rissues.htm
Box.test（fit2$resid,lag=6,type="Ljung-Box",fitdf=1）
fore=predict（fit2,10,newxreg=（length（incom）+1）:
（length（incom）+10））
#疏系数模型
#例5.8
w=ts（read.csv（"w.csv"）,start=1917,freq=1）
w=w[,1]
par（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（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.9
ue=ts（scan（"unemployment.txt"）,start=1962,f=4）
升级会员 