时间序列分析R语言程序.docx
《时间序列分析R语言程序.docx》由会员分享，可在线阅读，更多相关《时间序列分析R语言程序.docx（17页珍藏版）》请在冰豆网上搜索。
时间序列分析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次化学反应的过程数据进行拟合
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（iin1:
3）
{x[1]=x[1]
x[i+1]=0.25*x[i+1]+0.75*x[i]}
##错误于`[.data.frame`（x,i+1）:
undefinedcolumnsselected
##此外:
警告信息：
##InOps.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.7ARIMA
#例5.8疏系数模型妇女
Data1.15=read.csv（"C:
\\U
升级会员 