时间序列分析R语言程序

#例2.1绘制1964—— 1999年中国年纱产量序列时序 图(数据见附录1.2 )

Data1.2=read.csv(\附录

1.2.csv\

没有标题用F 如果有标题，用T;

#例 2.5对1950―― 1998年北京市城乡居民定期储蓄 所占比例序列的平稳性与纯随机性进行检验 Data1.5=read.csv(\

附录 1.5.csv\

plot(Data1.5,type='o',xlim=c(1950,2010),ylim=c (60,100)) 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(\附录

1.3.csv\

tdat1.3=as.vector(t(as.matrix(Data1.3)))[1:168 ]#矩阵转置转向量 plot(tdat1.3,type=T) #例 2.2 续 acf(tdat1.3) # 把字去掉

pacf(tdat1.3)

#例 2.3绘制1949―― 1998年北京市每年最高气温序 列时序图

Data1.4=read.csv(\附录

1.4.csv\

plot(Data1.4,type='o') ##不会定义坐标轴 #例 2.3 续

tdat1.4=Data1.4[,2] a1.4=acf(tdat1.4) #例 2.3 续

Box.test(tdat1.4,type=\Box.test(tdat1.4,type=\

#例 2.4随机产生1000个服从标准正态分布的白噪声 序列观察值，并绘制时序图 Data2.4=rnorm(1000,0,1) Data2.4

plot(Data2.4,type=T) #例 2.4 续 a2.4=acf(Data2.4) #例 2.4 续

Box.test(Data2.4,type=\Box.test(Data2.4,type=\

tdat1.5=Data1.5[,2] a1.5=acf(tdat1.5) #白噪声检验

Box.test(tdat1.5,type=\Box.test(tdat1.5,type=\

#例2.5续选择合适的ARMA模型拟合序列 acf(tdat1.5) pacf(tdat1.5)

#根据自相关系数图和偏自相关系数图可以判断为 (1)模型

#例2.5续P81 口径的求法在文档上 #P83

arima(tdat1.5,order=c(1,0,0),method=\似然估计

ar1= arima(tdat1.5,order=c(1,0,0),method=\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=\Box.test(ev,type=\#例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=T,col =1:2)

AR极大

lin es(U,col=\lin es(L,col=\#例3.1.1 例3.5 例3.5续

#方法一 plot.ts(arima.sim(n=100,list(ar=0.8))) #方法二 xO=ru nif(1) x=rep(0,1500) x[1]=0.8*x0+rnorm(1) for(i in 2:le ngth(x))

{x[i]=0.8*x[i-1]+rnorm(1)} plot(x[1:100],type=T) acf(x) pacf(x)

##拟合图没有画出来 #例 3.1.2 xO=ru nif(1) x=rep(0,1500) x[1]=-1.1*xO+rnorm(1) for(i in 2:le ngth(x))

{x[i]=-1.1*x[i-1]+rnorm(1)} plot(x[1:100],type=T) acf(x) pacf(x)

#例 3.1.3 方法一

plot.ts(arima.sim( n=100,list(ar=c(1,-0.5)))) #方法二 xO=ru nif(1) x1=ru nif(1) x=rep(0,1500) x[1]=x1

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

{x[i]=x[i-1]-0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=T) acf(x) pacf(x) #例 3.1.4 xO=ru nif(1) x1=ru nif(1) x=rep(0,1500) x[1]=x1

x[2]=x1+0.5*x0+rnorm(1) for(i in 3:le ngth(x))

{x[i]=x[i-1]+0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=T) acf(x) pacf(x) 又一个式子

xO=ru nif(1) x1=ru nif(1) x=rep(0,1500) x[1]=x1

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

{x[i]=-x[i-1]-0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=T) acf(x) pacf(x) #均值和方差 smu=mea n(x) svar=var(x)

#例3.2求平稳AR (1)模型的方差 例3.3 mu=0

mvar=1/(1-0.8A2) # 书上 51 页 #总体均值方差

cat(\mean and var

are\#样本均值方差

cat(\#例题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=T) acf(x) pacf(x) 法二

x=rep(0:1000) for(i in 1:1000)

{x[i]=rnorm[i]-2*rnorm[i-1]} plot(x,type=T) acf(x) pacf(x) #3.6.2 法一：

x=arima.sim (n=1000,list(ma=-0.5)) plot.ts(x,type=T) acf(x) pacf(x) 法二

x=rep(0:1000) for(i in 1:1000)

{x[i]=rnorm[i]-0.5*rnorm[i-1]} plot(x,type=T) 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=T) acf(x) pacf(x) 法二： x=rep(0:1000) for(i in 1:1000)

{x[i]=rnorm[i]-4/5*rnorm[i-1]+16/25*rnorm[i-2] } plot(x,type=T) 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=ru nif(1) x=rep(0,1000)

x[1]=0.5*x0+rnorm(1)-0.8*rnorm(1) for(i in 2:le ngth(x))

{x[i]=0.5*x[i-1]+rnorm(1)-0.8*rnorm(1)} plot(x,type=T) acf(x) pacf(x)

##图和书上不一样

#法二

x=arima.sim( n=1000,list(ar=0.5,ma=-0.8)) acf(x) pacf(x) #图和书上一样

#例3.8选择合适的 ARMA模型拟合加油站 57天的

OVERSHOR序列

Data1.6=read.csv(\

附录 1.6.csv\

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=\最 小二乘估计

ma1=arima(tdat1.6,order=c(0,0,1),method=\summary(ma1) ev=ma1$residuals acf(ev) pacf(ev)

##错误于 arima(tdat1.6, order = c(0, 0,

1),

method = \##'x'必需为数值

#例3.9选择合适的ARMAI型拟合1880—— 1985年全 球气温改变差值差分序列 ##没有数据

#例3.10 例3.11例3.12##矩估计

#例3.13对等时间间隔的连续 70次化学反应的过程数 据进行拟合

Data1.8=read.csv(\

附录 1.8.csv\

tdat1.8=as.vector(t(as.matrix(Data1.8)))[1:70] plot(tdat1.8,type='o')