时间序列分析R语言程序

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

Data1.2=read.csv(\附录1.2.csv\如果有标题，用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(\附录1.3.csv\

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(\附录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='l') #例2.4续

a2.4=acf(Data2.4) #例2.4续

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

#例2.5对1950——1998年北京市城乡居民定期储蓄所占比例序列的平稳性与纯随机性进行检验

Data1.5=read.csv(\附录1.5.csv\

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=\Box.test(tdat1.5,type=\#例2.5续选择合适的ARMA模型拟合序列 acf(tdat1.5) pacf(tdat1.5)

#根据自相关系数图和偏自相关系数图可以判断为AR（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=\lines(U,col=\lines(L,col=\

#例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(i in 2:length(x))

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

##拟合图没有画出来

#例3.1.2 x0=runif(1) x=rep(0,1500)

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

{x[i]=-1.1*x[i-1]+rnorm(1)} plot(x[1:100],type=\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(i in 3:length(x))

{x[i]=x[i-1]-0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=\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(i in 3:length(x))

{x[i]=x[i-1]+0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=\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(i in 3:length(x))

{x[i]=-x[i-1]-0.5*x[i-2]+rnorm(1)} plot(x[1:100],type=\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(\

#样本均值方差

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='l') acf(x) pacf(x) 法二

x=rep(0:1000) for(i in 1: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(i in 1: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(i in 1: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(i in 2: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(\附录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选择合适的ARMA模型拟合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')

#例3.14AR（2）例3.15AR（3）例3.16AR（3）模型的预测

#如果考得话就先。。。。。。

#例4.1线性拟合消费支出数据

Data4.1=read.csv(\例题4.1.csv\tdat4.1=Data4.1[,2] plot(Data4.1,type='o') t=1:40

lm4.1=lm(tdat4.1~t) #线性拟合

summary(lm4.1) #返回拟合参数的统计量 coef(lm4.1) #返回被估计的系数 fit4.1=fitted(lm4.1) #返回模拟值 residuals(lm4.1) #返回残差值 plot(tdat4.1,type='o') #画时序图

lines(fit4.1,col=\画拟合图

#例4.2 曲线拟合 上海证劵交易所

Data1.9=read.csv(\附录1.9.csv\

tdat1.9=as.vector(t(as.matrix(Data1.9)))[1:130]#矩阵转置转向量

plot(tdat1.9,type='l') t=1:130 t2=t^2

m1.9=lm(tdat1.9~t+t2) ## 一道矩阵就出毛病

#例4.3简单移动平均法 x4.3=c(5,5.4,5.8,6.2) x4.3

y4.3=filter(x4.3,rep(1/4,4),sides=1) y4.3

for(i in 1:3) {x[1]=x[1]

x[i+1]=0.25*x[i+1]+0.75*x[i]}

##错误于`[.data.frame`(x, i + 1) : undefined columns selected

##此外: 警告信息：

##In Ops.factor(left, right) : * 对因子没有意义

#例4.4指数平滑法 ##做不出来

#例4.5 ##略略

#例4.6季节效应分析

#例4.7综合分析 中国社会消费品零售总额序列 Data1.11=read.csv(\附录1.11.csv\第一行是标签，所以是T tdat1=as.matrix(Data1.11[,2:9]) # 横向全部读取，纵向读取2至9列

tdat1.11=as.vector(tdat1)

plot(1:length(tdat1.11),tdat1.11,type='o') #画时序图，先是横坐标，后是纵坐标

md=mean(tdat1.11)#求总的均值 md

seaind=apply(tdat1,1,'mean')/md #求季节因子 seaind

plot(seaind,type='b') # 季节指数图

noseandat=tdat1.11/seaind #消除季节因子的影响

plot(1:length(tdat1.11),noseandat,\消除季节因子之后的散点图

lindat=data.frame(x=1:length(noseandat),y=noseandat) m1=lm(y~x,data=lindat) #一元线性回归拟合 summary(m1)

t=1:96

that=1015.5222+20.9318*t

plot(1:length(tdat1.11),noseandat,'p')

lines(that,type='l') #拟合图和原来的图画在一起 #残差检验

ev=noseandat-that#计算残差 ev

plot(ev) #残差图 t=97:108

that=983.5601+21.5908*t q=that*seaind s=c(tdat1.11,q)

plot(1:108,s,type='b') abline(v=96)

#例5.1 差分运算

Data1.2=read.csv(\附录1.2.csv\如果有标题，用T；没有标题用F

x=Data1.2 plot(x,type='o') dx=diff(x[,2]) plot(dx,type='o')

#例5.2二阶差分 北京市民车辆拥有量序列

Data1.12=read.csv(\附录1.12.csv\x=Data1.12 plot(x,type='o')

dx=diff(x[,2]) #一届差分 plot(dx,type='b')

ddx=diff(x[,2],lag=1,difference=2) #二阶差分 plot(ddx,type='l')

#例5.2 又

Data1.12=read.csv(\附录1.12.csv\plot(Data1.12[,2],type='l') axis(1,at=c(1950,19999)) x=ts(Data1.12[,2])

dx=diff(x,lag=1,differences=1) #一阶差分 plot(Data1.12[-1,1],dx,type='o') d2x=diff(dx) #二阶差分

plot(Data1.12[-c(1,2),1],d2x,type='o') d2x=diff(x,differences=2) #二阶差分

plot(Data1.12[-c(1,2),1],d2x,type='o')

#例5.3 跳步差分 平均每头奶牛产奶量

Data1.13=read.csv(\附录1.13.csv\tdat1=as.matrix(Data1.13) tdat=t(tdat1) x=as.vector(tdat) x

plot(x,xaxt='n',type='o')

axis(1,at=seq(1,169,24),seq(1962,1976,2)) dx=diff(x) #一步差分 plot(dx,type='o')

d12x=diff(dx,lag=12) #12步差分 plot(d12x,type='o')

#例5.4 过差分

#例5.5 你和随机游走模型

r=rnorm(1000,sd=10) #以十位等差，在1：000之间随机抽取100个数据

xt=cumsum(r) #由随机游走公式的出的模型公式 plot(xt,type='l')# 随机游走的图形

dx=diff(xt) #做一阶差分

plot(dx,type='l') #一阶差分后的图形 m=mean(dx) #均值 sd=var(dx) #方差

Box.test(dx,lag=12,type='Ljung') #用统计量检验随机性

acf(dx) #自相关图

sj=arima(xt,order=c(0,1,0)) summary(sj)

#例5.5 你和随机游走模型 又 x=ts(cumsum(rnorm(1000,0,100))) ts.plot(x)

#例5.6 对中国农业实际国民收入指数进行建模ARIMA

Data1.14=read.csv(\附录1.14.csv\x=Data1.14

plot(x,type='o') #图5——10 dx=diff(x[,2]) plot(dx,type='o') acf(dx)

Box.test(dx,lag=6,type='Ljung-Box') Box.test(dx,lag=12,type='Ljung-Box') Box.test(dx,lag=18,type='Ljung-Box') pacf(dx)

m1=arima(x[,2],order=c(0,1,1),method=\m2=arima(dx,order=c(0,1,1),method=\ev1=m1$residuals ev2=m2$residuals plot(ev1,type='l') plot(ev2,type='l') acf(ev1) pacf(ev1) acf(ev2) pacf(ev2)

Box.test(ev1,lag=5,type='Ljung-Box') #检验残差的白噪声序列

Box.test(ev1,lag=11,type='Ljung-Box') Box.test(ev1,lag=17,type='Ljung-Box')

#例5.6续 做预测 ##没做好 px=predict(m1,n.ahead=10) plot(x,type='o',ylim=c(0,500)) lines(x[,1],x[,2]+1.96*sqrt(61.95)) lines(x[,1],x[,2]-1.96*sqrt(61.95)) ##图5——14没画出来

#例5.6续

m3=arima(x[,2],order=c(0,1,1),method=\

#例5.6续 p-171

plot(x,xlim=c(1950,1990),ylim=c(0,300),type='o') m1=lm(农业~年份,data=x) #变量为时间t的函数 summary(m1) ##？？？模型口径不会算 lines(x$年份,m1$fitted.value,col='red') #变量为一阶延迟 xt=x[,2] xy=xt[-1]

xx=xt[-length(xt)] m2=lm(xy~xx) summary(m2) m3=lm(xy~xx+0) summary(m3)

lines(x$年份[-1],m2$fitted.value,col='blue') #图5——29 #DW检验 library(lmtest)

dwtest(m2)#加载程序包 aa=dwtest(m1)

Dh=(1-aa$statistic/2)*sqrt(length(xt)-1)/(1-(length)-1)*0.009063 #Dh统计量 ev1=m1$residuals

plot(ev1,type='o')

m4=arima(ev1,order=c(2,0,0),fixed=c(NA,NA,0),trannsform.pars=F)

ev2=m3$residuals plot(ev2,type='o')

m5=arima(ev2,order=c(2,0,0),fixed=c(NA,0),trannsform.pars=F)

m6=arima(xt,order=c(0,1,1),xreg=1:length(xt),method='CSS')

#例5.7 ARIMA

#例5.8 疏系数模型 妇女

Data1.15=read.csv(\附录1.15.csv\x=Data1.15 plot(x,type='o') dx=diff(x[,2]) plot(dx,type='o') acf(dx) pacf(dx)

m=arima(x[,2],order=c(4,1,0),fixed=c(NA,0,0,NA),transform.pars=FALSE,method=\#不知道疏系数模型是怎么判断的 summary(m) ev=m$residuals

Box.test(ev,lag=6,type='Ljung-Box',fitdf=2) #阶数为原来的阶数减去参数的个数

Box.test(ev,lag=12,type='Ljung-Box',fitdf=2)

Box.test(ev,lag=18,type='Ljung-Box',fitdf=2) ##结果和书上不一样 #参数显著性检验 t1=0.2583/0.1159 2.228645

t2=0.3408/0.1225 2.782041 #t统计量

p1=pt(2.228645,df=57,lower.tail=F)*2

p2=pt(2.782041,df=57,lower.tail=F)*2 #p值不大对

#例5.9简单季节模型 德国工人失业率

Data1.16=read.csv(\附录1.16.csv\

tdat1.16=as.vector(t(as.matrix(Data1.16)))[1:120] #绘制时序图

x=ts(tdat1.16,start=1961,f=4) plot(x,type='o') #差分平稳化

dx1=diff(x)

plot(dx1,type='o') dx=diff(dx1,lag=4)

plot(dx,type='o',ylim=c(-2,2))

Box.test(dx,lag=6,type='Ljung-Box') Box.test(dx,lag=12,type='Ljung-Box')

Box.test(dx,lag=18,type='Ljung-Box') #表5——6 差分序列具有很强的相关信息 #模型拟合 acf(dx) pacf(dx)

m=arima(x,order=c(4,1,0),fixed=c(NA,0,0,NA),transform.pars=FALSE,include.mean=F,method=\coef(m) ##模型拟合的不对 #参数估计与检验 ##不会 ev1=m$residuals

Box.test(ev,lag=4,type='Ljung-Box') Box.test(ev,lag=10,type='Ljung-Box')

#例5.10乘积季节模型

Data1.17=read.csv(\附录1.17.csv\

tdat1.17=as.vector(t(as.matrix(Data1.17)))[1:408] x=ts(tdat1.17,start=1948,f=12) plot(x,type='l',ylim=c(0,4000)) #差分平稳化 dx1=diff(x)

plot(dx1,type='l',ylim=c(-400,500)) dx=diff(dx1,lag=12)

plot(dx,type='l',ylim=c(-400,500)) #图5——25 acf(dx) pacf(dx)

#例5.11 异方差的性质

Data1.18=read.csv(\附录1.18.csv\Data1.18

tdat1.18=as.vector(t(as.matrix(Data1.18)))[1:100] x=ts(tdat1.18,start=1962,f=12) plot(x,type='l',ylim=c(0,0.0065)) dx1=diff(x)

plot(dx1,type='l',ylim=c(-0.0010,0.0012)) #图5-38 acf(dx1) pacf(dx1)

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库时间序列分析R语言程序在线全文阅读。