应⽤时间序列分析第三、四章代码
#unit3
library("timeSeries")
#例3-1
x1<-arima.sim(n=100,list(ar=0.8))
x3<-arima.sim(n=100,list(ar=c(1,-0.5)))
e<-rnorm(100)
x2<-filter(e,filter =-1.1,method ="recursive")
x4<-filter(e,filter = c(1,0.5),method ="recursive")
ts.plot(x1)
ts.plot(x2)
ts.plot(x3)
ts.plot(x4)
#例3-5
x1<-arima.sim(n=1000,list(ar=0.8))
x2<-arima.sim(n=1000,list(ar=-0.8))
x3<-arima.sim(n=1000,list(ar=c(1,-0.5)))
x4<-arima.sim(n=1000,list(ar=c(-1,-0.5)))
acf(x1)
acf(x2)
acf(x3)
acf(x4)
#例3-5 续
#par(mfcol=c(2,2))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)
#例3-6
#par(mfcol=c(2,2))
x1<-arima.sim(n=1000,list(ma=-2))
运营的英文x2<-arima.sim(n=1000,list(ma=-0.5))
x3<-arima.sim(n=1000,list(ma=c(-4/5,16/25)))
x4<-arima.sim(n=1000,list(ma=c(-5/4,25/16)))
acf(x1)
acf(x2)
acf(x3)
张姓女孩名acf(x4)
#例3-6 续(2)
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)
束缚什么意思#例3-8
x<-arima.sim(n=1000,list(ar=0.5,ma=-0.8))
acf(x)
pacf(x)
#例3-9 1-8
#读⼊数据，并绘制时序图
#读⼊附录1-8一伏是多少天
#a<-read.table("E:/R/data/file8.csv",p=",",header = T)
x<-ts(X3_9$kilometer,start =1950)
plot(x)
#⽩噪声检验
for(i in1:2) st(x,type ="Ljung-Box",lag=6*i))
#绘制⾃相关图和偏⾃相关图
pacf(x)
#例3-10
#读⼊数据，并绘制时序图附录1-9
#overshort<-read.table("E:/R/data/file9.csv",p=",",header = T) overshort<-ts(X3_10$overshort,start =1)
#fulu1_9 timeries name
plot(overshort)
#⽩噪声检验
for(i in1:2) st(overshort,type ="Ljung-Box",lag=6*i)) #绘制⾃相关图和偏⾃相关图
acf(overshort)
什么头什么尾pacf(overshort)
#例3-11
#读⼊数据，并绘制时序图附录10
#b<-read.table("E:/R/data/file10.csv",p=",",header = T)
dif_x<-ts(diff(X3_11$change_temp),start =1880)
plot(dif_x)
#⽩噪声检验
for(i in1:2) st(dif_x,type ="Ljung-Box",lag=6*i))
#绘制⾃相关图和偏⾃相关图
acf(dif_x)
pacf(dif_x)
久开头的成语install.packages("zoo")
install.packages("forecast")
library(zoo)
library(forecast)
#例3-9系统⾃动定阶
auto.arima(x)
#例3-10系统⾃动定阶
auto.arima(overshort)
#例3-11系统⾃动定阶
auto.arima(dif_x)
# 参数估计
#例3-9续(1)附录8
#a<-read.table("E:/R/data/file8.csv",p=",",header = T) kilometer<-ts(X3_9$kilometer,start =1950)
kilometer.fit<-arima(kilometer,order = c(2,0,0),method ="ML") kilometer.fit
#例3-10续(1)
#overshort<-read.table("E:/R/data/file9.csv",p=",",header = T) #overshort<-ts(overshort)
overshort.fit<-arima(overshort,order = c(0,0,1))
overshort.fit
#例3-11续(1)
#b<-read.table("E:/R/data/file10.csv",p=",",header = T)越南游记
dif_x<-ts(diff(X3_11$change_temp),start =1880)
dif_x.fit<-arima(dif_x,order = c(1,0,1))
dif_x.fit
#模型显著性检验
x<-ts(a$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0),method ="ML")
for(i in1:2) st(x.fit$residual,lag=6*i))
#例3-10续(2)
#overshort<-read.table("E:/R/data/file9.csv",p=",",header = T) overshort<-ts(overshort)
overshort.fit<-arima(overshort,order = c(0,0,1))
for(i in1:2) st(overshort.fit$residual,lag=6*i))
#例3-11续(2)
#b<-read.table("E:/R/data/file10.csv",p=",",header = T)
dif_x<-ts(diff(b$chang_temp),start =1880)
dif_x.fit<-arima(dif_x,order = c(1,0,1),method ="CSS")
for(i in1:2) st(dif_x.fit$residual,lag=6*i))
#参数显著性检验
#例3-9续(3)
#a<-read.table("E:/R/data/file8.csv",p=",",header = T)
x<-ts(x3_9$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0),method ="ML")
x.fit
#ar1系数显著性检验
t1<-0.7185/0.1083
pt(t1,df=56,lower.tail = F)
#ar2系数显著性检验
t2<-0.5294/0.1067
pt(t2,df=56,lower.tail = T)
徒子徒孙#ar3系数显著性检验
t0=11.0223/3.0906
pt(t0,df=56,lower.tail = F)
#模型优化
#例3-15
#读⼊数据，绘制时序图
#x<-read.table(file = "E:/R/data/file11.csv",p=",",header = T) y<-ts(X3_15$yield)
plot(y)
#序列⽩噪声检验
for(i in1:2) st(y,lag=6*i))
#绘制⾃相关图和偏⾃相关图
which(is.na(y))#查找缺失值
y[24]<- mean( = T)#缺失值处理
y[64]<- mean( = T)
acf(y)
pacf(y)
#拟合MA(2)模型
y.fit1<-arima(y,order = c(0,0,2))
y.fit1
#MA(2)模型显著性检验
for(i in1:2) st(y.fit1$residual,lag=6*i))#残差检验
#拟合AR(1)模型
y.fit2<-arima(y,order = c(1,0,0))
y.fit2
#AR(1)模型显著性检验
for(i in1:2) st(y.fit2$residual,lag=6*i))
#序列预测
library(forecast)
x<-ts(X3_9$kilometer,start=1950)
x.fit<-arima(x,order = c(2,0,0))#由于x的偏⾃相关明显2阶结尾，则选⽤AR(2) x.fore<-forecast(x.fit,h=5)#h：预测期数
x.fore
#系统默认输出预测图
plot(x.fore)
#个性化输出预测图
L1<-x.fore$fitted-1.96*sqrt(x.fit$sigma2)#fitted代表拟合值,
U1<-x.fore$fitted+1.96*sqrt(x.fit$sigma2)
L2<-ts(x.fore$lower[,2],start =2009)#预测值95%的置信区间下限
U2<-ts(x.fore$upper[,2],start =2009)#预测值95%的置信区间上限
c1<-min(x,L1,L2)
c2<-max(x,L2,U2)
plot(x,type ="p",pch=8,xlim = c(1950,2013),ylim = c(c1,c2))
lines(x.fore$fitted,col=2,lwd=2)#拟合曲线
lines(x.fore$mean,col=2,lwd=2)#添加预测曲线
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(U2,col=4,lty=2)
