第一讲 应用实例
? R的基本界面是一个交互式命令窗口,命令提示符是一个大于号,命令的结果马上
显示在命令下面。
? S命令主要有两种形式:表达式或赋值运算(用’<-’或者’=’表示)。在命令提示符后
键入一个表达式表示计算此表达式并显示结果。赋值运算把赋值号右边的值计算出来赋给左边的变量。
? 可以用向上光标键来找回以前运行的命令再次运行或修改后再运行。 ? S是区分大小写的,所以x和X是不同的名字。
我们用一些例子来看R软件的特点。假设我们已经进入了R的交互式窗口。如果没有打开的图形窗口,在R中,用:> x11() 可以打开一个作图窗口。然后,输入以下语句: x1 = 0:100
x2 = x1*2*pi/100 y = sin(x2)
plot(x2,y,type=\
这些语句可以绘制正弦曲线图。其中,“=”是赋值运算符。0:100表示一个从0到100 的等差数列向量。第二个语句可以看出,我们可以对向量直接进行四则运算,计算得到的x2 是向量x1的所有元素乘以常数2*pi/100的结果。从第三个语句可看到函数可以以向量为输入,并可以输出一个向量,结果向量y的每一个分量是自变量x2的每一个分量的正弦函数值。
plot(x2,y, type=\画图练习\好好练\轴\轴')
有关作图命令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(\ 注意字符串中的反斜杠。
例:计算6, 4, 7, 8,10的均值和标准差,把若干行命令保存在一个文本文件(比如C:\\1.R)中,然后用source函数来运行整个文件。 a<- c(10, 6, 4, 7, 8) b<-mean(a) c<-sd(a)
source(\
时间序列数据的输入 使用函数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=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 = \ #保存协方差 图的保存,单击选中图,在菜单栏选中“文件”,再选“另存为”。 同时显示多个图:用x11()命令生成一个空白图,再输入作图命令。
2. 同时绘制两组数据的时序图 d=read.csv(\double=ts(d,start=1964,freq=1)
plot(double, plot.type = \ #两组数据两个图 plot(double, plot.type = \ #两组数据一个图
plot(double, plot.type = \设置每组数据图的颜色、曲线类型)
3.产生服从正态分布的随机观察值
例题2.4 随机产生1000白噪声序列观察值 d=rnorm(1000,0,1) #个数1000 均值0 方差1 plot.ts(d)
4.纯随机性检验 例题2.3续
d=scan(\
temp=ts(d,freq=1,start=c(1949))
Box.test(temp, type=\
5.差分计算 x=1:10 y=diff(x) k步差分
?kxt?xt?xt?k 加入参数 lag=k
如计算x的3步差分为 y=diff(x, lag = 3)
pp?1p?1?x??x??xt?1加入参数differences = p ttp阶差分
2?xt??xt??xt?1
如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(\ #导入数据
prop=ts(d,start=1950,freq=1) #转化为时间序列数据 plot(prop) #作时序图
acf(prop,12) #作自相关图,拖尾
pacf(prop,12) #作偏自相关图,1阶截尾 Box.test(prop, type=\ #纯随机性检验,p值小于5%,序列为非白噪声 Box.test(prop, type=\arima(prop, order = c(1,0,0),method=\
#用AR(1)模型拟合,如参数method=\,估计方法为条件最小二乘法,用条件最小二乘法时,不显示AIC。
arima(prop, order = c(1,0,0),method=\用AR(1)模型拟合,不含截距项。 tsdiag(arima(prop, order = c(1,0,0),method=\ #对估计进行诊断,判断残差是否为白噪声
summary(arima(prop, order = c(1,0,0),method=\a=arima(prop, order = c(1,0,0),method=\r=a$residuals#用r来保存残差
Box.test(r,type=\对残差进行纯随机性检验
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=\
lines(L, col=\在时序图中作出95%置信区间
例题3.9
d=scan(\x=diff(d)
arima(x, order = c(1,0,1),method=\
tsdiag(arima(x, order = c(1,0,1),method=\
第一点:
对于第三讲中的例2.5,运行命令arima(prop, order = c(1,0,0),method=\之后,显示: Call:
arima(x = prop, order = c(1, 0, 0), method = \Coefficients:
ar1 intercept 0.6914 81.5509 s.e. 0.0989 1.7453
sigma^2 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.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值,如果题目没有特别要求,一般不需要手动计算。
第三点:
修正第三讲中的错误: