时间序列分析R语言程序.docx
- 文档编号:474270
- 上传时间:2022-10-10
- 格式:DOCX
- 页数:17
- 大小:19.84KB
时间序列分析R语言程序.docx
《时间序列分析R语言程序.docx》由会员分享,可在线阅读,更多相关《时间序列分析R语言程序.docx(17页珍藏版)》请在冰豆网上搜索。
时间序列分析R语言程序
时间序列分析R语言程序
#例2.1绘制1964——1999年中国年纱产量序列时序图(数据见附录1.2)
Data1.2=read.csv("C:
\\Users\\Administrator\\Desktop\\附录1.2.csv",header=T)#如果有标题,用T;没有标题用F
plot(Data1.2,type='o')
#例2.1续
tdat1.2=Data1.2[,2]
a1.2=acf(tdat1.2)
#例2.2绘制1962年1月至1975年12月平均每头奶牛产奶量序列时序图(数据见附录1.3)
Data1.3=read.csv("C:
\\Users\\Administrator\\Desktop\\附录1.3.csv",header=F)
tdat1.3=as.vector(t(as.matrix(Data1.3)))[1:
168]#矩阵转置转向量
plot(tdat1.3,type='l')
#例2.2续
acf(tdat1.3)#把字去掉
pacf(tdat1.3)
#例2.3绘制1949——1998年北京市每年最高气温序列时序图
Data1.4=read.csv("C:
\\Users\\Administrator\\Desktop\\附录1.4.csv",header=T)
plot(Data1.4,type='o')
##不会定义坐标轴
#例2.3续
tdat1.4=Data1.4[,2]
a1.4=acf(tdat1.4)
#例2.3续
Box.test(tdat1.4,type="Ljung-Box",lag=6)
Box.test(tdat1.4,type="Ljung-Box",lag=12)
#例2.4随机产生1000个服从标准正态分布的白噪声序列观察值,并绘制时序图
Data2.4=rnorm(1000,0,1)
Data2.4
plot(Data2.4,type='l')
#例2.4续
a2.4=acf(Data2.4)
#例2.4续
Box.test(Data2.4,type="Ljung-Box",lag=6)
Box.test(Data2.4,type="Ljung-Box",lag=12)
#例2.5对1950——1998年北京市城乡居民定期储蓄所占比例序列的平稳性与纯随机性进行检验
Data1.5=read.csv("C:
\\Users\\Administrator\\Desktop\\附录1.5.csv",header=T)
plot(Data1.5,type='o',xlim=c(1950,2010),ylim=c(60,100))
tdat1.5=Data1.5[,2]
a1.5=acf(tdat1.5)
#白噪声检验
Box.test(tdat1.5,type="Ljung-Box",lag=6)
Box.test(tdat1.5,type="Ljung-Box",lag=12)
#例2.5续选择合适的ARMA模型拟合序列
acf(tdat1.5)
pacf(tdat1.5)
#根据自相关系数图和偏自相关系数图可以判断为AR
(1)模型
#例2.5续P81口径的求法在文档上
#P83
arima(tdat1.5,order=c(1,0,0),method="ML")#极大似然估计
ar1=arima(tdat1.5,order=c(1,0,0),method="ML")
summary(ar1)
ev=ar1$residuals
acf(ev)
pacf(ev)
#参数的显著性检验
t1=0.6914/0.0989
p1=pt(t1,df=48,lower.tail=F)*2
#ar1的显著性检验
t2=81.5509/1.7453
p2=pt(t2,df=48,lower.tail=F)*2
#残差白噪声检验
Box.test(ev,type="Ljung-Box",lag=6,fitdf=1)
Box.test(ev,type="Ljung-Box",lag=12,fitdf=1)
#例2.5续P94预测及置信区间
predict(arima(tdat1.5,order=c(1,0,0)),n.ahead=5)
tdat1.5.fore=predict(arima(tdat1.5,order=c(1,0,0)),n.ahead=5)
U=tdat1.5.fore$pred+1.96*tdat1.5.fore$se
L=tdat1.5.fore$pred-1.96*tdat1.5.fore$se
plot(c(tdat1.5,tdat1.5.fore$pred),type="l",col=1:
2)
lines(U,col="blue",lty="dashed")
lines(L,col="blue",lty="dashed")
#例3.1.1例3.5例3.5续
#方法一plot.ts(arima.sim(n=100,list(ar=0.8)))
#方法二
x0=runif
(1)
x=rep(0,1500)
x[1]=0.8*x0+rnorm
(1)
for(iin2:
length(x))
{x[i]=0.8*x[i-1]+rnorm
(1)}
plot(x[1:
100],type="l")
acf(x)
pacf(x)
##拟合图没有画出来
#例3.1.2
x0=runif
(1)
x=rep(0,1500)
x[1]=-1.1*x0+rnorm
(1)
for(iin2:
length(x))
{x[i]=-1.1*x[i-1]+rnorm
(1)}
plot(x[1:
100],type="l")
acf(x)
pacf(x)
#例3.1.3
方法一
plot.ts(arima.sim(n=100,list(ar=c(1,-0.5))))
#方法二
x0=runif
(1)
x1=runif
(1)
x=rep(0,1500)
x[1]=x1
x[2]=x1-0.5*x0+rnorm
(1)
for(iin3:
length(x))
{x[i]=x[i-1]-0.5*x[i-2]+rnorm
(1)}
plot(x[1:
100],type="l")
acf(x)
pacf(x)
#例3.1.4
x0=runif
(1)
x1=runif
(1)
x=rep(0,1500)
x[1]=x1
x[2]=x1+0.5*x0+rnorm
(1)
for(iin3:
length(x))
{x[i]=x[i-1]+0.5*x[i-2]+rnorm
(1)}
plot(x[1:
100],type="l")
acf(x)
pacf(x)
又一个式子
x0=runif
(1)
x1=runif
(1)
x=rep(0,1500)
x[1]=x1
x[2]=-x1-0.5*x0+rnorm
(1)
for(iin3:
length(x))
{x[i]=-x[i-1]-0.5*x[i-2]+rnorm
(1)}
plot(x[1:
100],type="l")
acf(x)
pacf(x)
#均值和方差
smu=mean(x)
svar=var(x)
#例3.2求平稳AR
(1)模型的方差例3.3
mu=0
mvar=1/(1-0.8^2)#书上51页
#总体均值方差
cat("populationmeanandvarare",c(mu,mvar),"\n")
#样本均值方差
cat("samplemeanandvarare",c(mu,mvar),"\n")
#例题3.4
svar=(1+0.5)/((1-0.5)*(1-1-0.5)*(1+1-0.5))
#例题3.6MA模型自相关系数图截尾和偏自相关系数图拖尾
#3.6.1
法一:
x=arima.sim(n=1000,list(ma=-2))
plot.ts(x,type='l')
acf(x)
pacf(x)
法二
x=rep(0:
1000)
for(iin1:
1000)
{x[i]=rnorm[i]-2*rnorm[i-1]}
plot(x,type='l')
acf(x)
pacf(x)
#3.6.2
法一:
x=arima.sim(n=1000,list(ma=-0.5))
plot.ts(x,type='l')
acf(x)
pacf(x)
法二
x=rep(0:
1000)
for(iin1:
1000)
{x[i]=rnorm[i]-0.5*rnorm[i-1]}
plot(x,type='l')
acf(x)
pacf(x)
##错误于rnorm[i]:
类别为'closure'的对象不可以取子集
#3.6.3
法一:
x=arima.sim(n=1000,list(ma=c(-4/5,16/25)))
plot.ts(x,type='l')
acf(x)
pacf(x)
法二:
x=rep(0:
1000)
for(iin1:
1000)
{x[i]=rnorm[i]-4/5*rnorm[i-1]+16/25*rnorm[i-2]}
plot(x,type='l')
acf(x)
pacf(x)
##错误于x[i]=rnorm[i]-4/5*rnorm[i-1]+16/25*rnorm[i-2]:
##更换参数长度为零
#例3.6续根据书上64页来判断
#例3.7拟合ARMA(1,1)模型,x(t)-0.5x(t-1)=u(t)-0.8*(u-1),并直观观察该模型自相关系数和偏自相关系数的拖尾性。
#法一:
x0=runif
(1)
x=rep(0,1000)
x[1]=0.5*x0+rnorm
(1)-0.8*rnorm
(1)
for(iin2:
length(x))
{x[i]=0.5*x[i-1]+rnorm
(1)-0.8*rnorm
(1)}
plot(x,type='l')
acf(x)
pacf(x)
##图和书上不一样
#法二
x=arima.sim(n=1000,list(ar=0.5,ma=-0.8))
acf(x)
pacf(x)
#图和书上一样
#例3.8选择合适的ARMA模型拟合加油站57天的OVERSHORT序列
Data1.6=read.csv("C:
\\Users\\Administrator\\Desktop\\附录1.6.csv",header=F)
tdat1.6=as.vector(t(as.matrix(Data1.6)))[1:
57]
plot(tdat1.6,type='o')
acf(tdat1.6)
pacf(tdat1.6)#把字去掉
arima(tdat1.6,order=c(0,0,1),method="CSS")#最小二乘估计
ma1=arima(tdat1.6,order=c(0,0,1),method="CSS")
summary(ma1)
ev=ma1$residuals
acf(ev)
pacf(ev)
##错误于arima(tdat1.6,order=c(0,0,1),method="CSS"):
##'x'必需为数值
#例3.9选择合适的ARMA模型拟合1880——1985年全球气温改变差值差分序列
##没有数据
#例3.10例3.11例3.12##矩估计
#例3.13对等时间间隔的连续70
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 时间 序列 分析 语言 程序