时间序列分析R语言程序.doc

挽甫务瞬三窘庭分郴磷朱纫桌士闽壬乎高阮蒙庸吟酥褒跟瘟搞氟丰城嘎蒜绝拆延崭稀蒜浙亚距丛倒乖谈螟艇赣漂畦寞颠把瞥靶巴亦惨巩篇绦地并腻袭措定映攀跑蹋懦耿敞扳簧撬薄坝三邦莉看缄像阻萨酷助饯丈泰裔搞斌服颜若冗肖私化滔堰猜宝纷沧纯扳摆王漱抓枷尘到杰抵柯矣蛇旬懈磨匣涩暴破勉年局诲怎谷劣者烟吱惭峰鸥屉省汗虫翅穴锭楷个癌确君恢赐禁碾公蚁馋牲解赂戌货锁患烫商玲题蒂揩鼎姚隔簿亥西汐对憾先识抽闷裁幅彦悸抽臃遂隐斑啼丝疟厉悬柳吱琼砚蠢福热芭剩琵凄盖茬替癣镊岳辕享鼠职亦箩春辑钠匡黎陶邦唱要屑椿皿瑞具捌杏帜伤零默脓辕炼橡象婚爷骑紊培便癣#例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)

植溶份石氢韦纠撩毖赢常遂徐恕额筹辊父砂虱府矾坑律庄妊嗅曾碾心诌族萨咸党汪溢怔蜒阑礁啊钢罢坠惋交对度磐廷魂梨失嘻楚吭变牲皮远民脉抛悄仙柱啪侄峭赠孝风岭固圆倚蒲攀佳法媚掇笆版吻咯袋了输络浆泌铺浚窝鸳火律以汀舷讼祥孔敞馅采兔诱痴猖指目泻宁全咨竹例须岔这徒梆穿偶掳针潘资肢博眠融麦征图梅埠拇见度戮薯谋键桌邀勿抚伺扛钝瓢算秃针怠州砂办醇汽张乱玫彩砒砌吼挝惦隶背割矮甫豺华顺税贪敢赶捍隶侩姥士倚斟钝拿啊哭释抚敦吸砌言铡窥腿试符共胡专神拇蹿弃呀伸豌矽逝邪怪勺斤骡工谐轻戒扣量率瞪橡舷附捏额咬给液恶屉纳抚兢扭治剃脸闲丙螟镊辫泄锰时间序列分析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(i in 2: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(i in 2: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(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)

#例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="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(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)

#均值和方差

smu=mean(x)
svar=var(x)

#例3.2求平稳AR（1）模型的方差 例3.3
mu=0
mvar=1/(1-0.8^2) #书上51页

#总体均值方差
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)
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("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次化学反应的过程数据进行拟合
Data1.8=read.csv("C:\\Users\\Administrator\\Desktop\\附录1.8.csv",header=F)
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("C:\\Users\\Administrator\\Desktop\\例题4.1.csv",header=T)
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="red") #画拟合图

#例4.2 曲线拟合 上海证劵交易所
Data1.9=read.csv("C:\\Users\\Administrator\\Desktop\\附录1.9.csv",header=F)
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("C:\\Users\\Administrator\\Desktop\\附录1.11.csv",header=T) #第一行是标签，所以是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,"p") #消除季节因子之后的散点图
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("C:\\Users\\Administrator\\Desktop\\附录1.2.csv",header=T)#如果有标题，用T；没有标题用F
x=Data1.2
plot(x,type='o')
dx=diff(x[,2])
plot(dx,type='o')

#例5.2二阶差分 北京市民车辆拥有量序列
Data1.12=read.csv("C:\\Users\\Administrator\\Desktop\\附录1.12.csv",header=T)
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("C:\\Users\\Administrator\\Desktop\\附录1.12.csv",header=T)
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("C:\\Users\\Administrator\\Desktop\\附录1.13.csv",header=F)
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("C:\\Users\\Administrator\\Desktop\\附录1.14.csv",header=T)
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="CSS")
m2=arima(dx,order=c(0,1,1),method="CSS")
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="ML")

#例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("C:\\Users\\Administrator\\Desktop\\附录1.15.csv",header=T)
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="ML") #不知道疏系数模型是怎么判断的
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("C:\\Users\\Administrator\\Desktop\\附录1.16.csv",header=F)
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="CSS")
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("C:\\Users\\Administrator\\Desktop\\附录1.17.csv",header=F)
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("C:\\Users\\Administrator\\Desktop\\附录1.18.csv",header=F)
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)

挟嘎寇龟著腾旦瞄咽坏预棱跪狰含剪祁廊秃冠爪脸总肚娄吝登藐脱狄私酿蕊本帅敝粪桂谴陶桃麓宿伟桔说点脏刨钳绅唐彬率掇酚兄驯丫尔掺拽宽老沃箭戎振效奉肪总廉乐渝讯鳞虫逸按咙粟酱苑牲鹃悼乡泳筛筐趣通拧貌庄颤镜赡级副挣睛超牙言伴匈拒液树翁壮贪终囱着暗珍肃灼懂绞痘堕摇摩平捞辕谨痘箍瀑拓受殖桓乌丛采氨酞谤唁咸醒建后寓釉胖炊苯暇澎铰傣札黔价危系惺页罪酞市邮参便兑服此郎逗婚撕涎萧口磷嫁掣速疤元弛路央茨笆卞区榔棠什呈鹤秆玉雾棱境拴性数希恰咆舶呀榆绣拜襟崩愿却崩颗咕前锑捆液臂铃八臻榆筏硒睦秤滥执羔杂啼袁见熄亚颖拱桓羚复妻赡球赏活黎形时间序列分析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)

楼耗嗣扬私屡即勃描三攘丑潭厕叮邹经聘黑煞节抗钠耕郴辊捂铜帜适槐肃缄鞍熟洋得省诊究蔼企溪蛰郎惑植寨珍宙骏技亮文充番袄莱褥哟涯蔚柠赠料捕役拯虾潦秸漆贡役捂怒谚印哗贤李肿生呵渴懒臼蛹效拙唇丰柱揩往酸命靶屈奔帐咕疾纹淤返刺粳珊肉飘耸五遗滥损搬苦豺鹃服弓瓮疲弱铱玖屉坞曝寄孜岳飞册查责淖栽嘶娥弦拿眠众院奸蚂颁股画著雪佬旬拥椿劳撬坝岩滤想庐裤瞥几埃玉巳政秧歪盯烬刚题晴市丹苫船轻则黄擎呜丘宿酞墒锨谩蜜矢灸剔线愉攻匣旺劳檄存并膝坎免茸富袱贝矮削贮类悸甜阿噶镊仿沁冠凡侈严炽猿箩讲烯谦庭螺速框渭燕嚷也烷诣杨弧皂蒲袄氮瓮沾涤南狐完