时间序列分析R语言程序

#例 2.1 绘制 1964—— 1999 年中国年纱产量序列时序

图（数据见附录 1.2）

附录 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）

附录 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 年北京市每年最高气温序

列时序图

附录 1.4.csv",header=T)

plot(Data1.4,type='o')附录 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

#例 2.3 续p2=pt(t2,df=48,lower.tail=F)*2

tdat1.4=Data1.4[,2]#残差白噪声检验

a1.4=acf(tdat1.4)Box.test(ev,type="Ljung-Box",lag=6,fitdf=1)

#例 2.3 续Box.test(ev,type="Ljung-Box",lag=12,fitdf=1)

Box.test(tdat1.4,type="Ljung-Box",lag=6)#例 2.5 续 P94 预测及置信区间

Box.test(tdat1.4,type="Ljung-Box",lag=12)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.ahea #例 2.4 随机产生1000 个服从标准正态分布的白噪声d=5)

序列观察值，并绘制时序图

Data2.4=rnorm(1000,0,1)U=tdat1.5.fore$pred+1.96*tdat1.5.fore$se

Data2.4L=tdat1.5.fore$pred-1.96*tdat1.5.fore$se

plot(Data2.4,type='l')

#例 2.4 续plot(c(tdat1.5,tdat1.5.fore$pred),type="l",col=1:2)

a2.4=acf(Data2.4)lines(U,col="blue",lty="dashed")

#例 2.4 续lines(L,col="blue",lty="dashed")

Box.test(Data2.4,type="Ljung-Box",lag=6)

Box.test(Data2.4,type="Ljung-Box",lag=12)#例 3.1.1 例 3.5例3.5续

#方法一 plot.ts(arima.sim(n=100,list(ar=0.8)))

#例 2.5 对 1950—— 1998 年北京市城乡居民定期储蓄#方法二

所占比例序列的平稳性与纯随机性进行检验x0=runif(1)

x=rep(0,1500)

x[1]=0.8*x0+rnorm(1)x[1]=x1

for(i in 2:length(x))x[2]=-x1-0.5*x0+rnorm(1)

{x[i]=0.8*x[i-1]+rnorm(1)}for(i in 3:length(x))

plot(x[1:100],type="l"){x[i]=-x[i-1]-0.5*x[i-2]+rnorm(1)}

acf(x)plot(x[1:100],type="l")

pacf(x)acf(x)

##拟合图没有画出来pacf(x)

#例 3.1.2#均值和方差

x0=runif(1)

x=rep(0,1500)smu=mean(x)

x[1]=-1.1*x0+rnorm(1)svar=var(x)

for(i in 2:length(x))

{x[i]=-1.1*x[i-1]+rnorm(1)}#例 3.2 求平稳 AR （ 1）模型的方差例 3.3 plot(x[1:100],type="l")mu=0

acf(x)mvar=1/(1-0.8^2) # 书上 51 页

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(i in 3:length(x))

{x[i]=x[i-1]-0.5*x[i-2]+rnorm(1)}

plot(x[1:100],type="l")

acf(x)

pacf(x)#总体均值方差

cat("population mean and var are",c(mu,mvar),"\n")

#样本均值方差

cat("sample mean and var are",c(mu,mvar),"\n")

#例题 3.4

svar=(1+0.5)/((1-0.5)*(1-1-0.5)*(1+1-0.5))

#例题 3.6 MA 模型自相关系数图截尾和偏自相关系数图拖尾

#3.6.1

法一：

x=arima.sim(n=1000,list(ma=-2))

plot.ts(x,type='l')

acf(x)

#例 3.1.4pacf(x)

x0=runif(1)

x1=runif(1)法二

x=rep(0,1500)x=rep(0:1000)

x[1]=x1for(i in 1:1000)

x[2]=x1+0.5*x0+rnorm(1){x[i]=rnorm[i]-2*rnorm[i-1]}

for(i in 3:length(x))plot(x,type='l')

{x[i]=x[i-1]+0.5*x[i-2]+rnorm(1)}acf(x)

plot(x[1:100],type="l")pacf(x)

acf(x)

pacf(x)#3.6.2

法一：

又一个式子x=arima.sim(n=1000,list(ma=-0.5)) x0=runif(1)plot.ts(x,type='l')

x1=runif(1)acf(x)

x=rep(0,1500)pacf(x)