第三次上机实验：模型定阶与参数估计
作业内容参见：http://wenku.baidu.com/view/c128c60702020740be1e9b14.html
参考解答：
# 3、根据之前上机实验课的知识模拟产生下列时间序列的10000个观测数据，其中，构成白噪声序列的随机变量服从标准正态分布
set.seed(100)
ar1.sim<-arima.sim(list(order=c(1,0,0), ar=-0.65), n=10000, rand.gen=rnorm)
ar2.sim<-arima.sim(list(order=c(2,0,0), ar=c(0.5,0.3)), n=10000, rand.gen=rnorm)
ma3.sim<-arima.sim(list(order=c(0,0,1), ma=-0.5), n=10000, rand.gen=rnorm)
ma4.sim<-arima.sim(list(order=c(0,0,2), ma=c(-0.65,-0.24)), n=10000, rand.gen=rnorm)
arma5.sim<-arima.sim(list(order=c(1,0,1), ar=0.8,ma=-0.4), n=10000, rand.gen=rnorm)
#3、 序列（1）的真实阶数为1，分别用AR(1)模型和AR(2)模型对观测数据进行拟合（当然，之前我们应该首先观察数据图形粗略判断平稳性，观察数据的样本ACF及样本PACF对模型形式有初步判断。由于此处我们已知该时间序列数据的生成模型，因此将这一步骤省略，但是实际分析问题时是不可省的），观察其AIC值。即输入下列语句：
x<-ar1.sim
#平稳性检验
par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture1
library(tseries)
adf.test(x)#p-value = 0.01,stationary
#数据拟合
result1<-arima(x,order=c(1,0,0))
result2<-arima(x,order=c(2,0,0))
result11<-arima(x,order=c(1,0,0),include.mean=F)
result21<-arima(x,order=c(2,0,0),include.mean=F)
#4、 记录上述两个模型的AIC值，确定拟合模型的阶数并与模型实际阶数相比较。
AIC(result1)#28286.07
AIC(result2)#28285.69
AIC(result11)#28284.25
AIC(result21)#28283.88
#单从aic的值来判断，第一个模型用ar（2）模拟似乎更好，这一点与模型并不符合，包括与pacf图像的观察结果来看也不符合。但是两个模型的aic的值差异并不大。所以单从aic来定阶也未必那么科学。我们不妨来看看ar（2）模型是否存在过拟合
r<-result1$residuals
plot(r,type='o')
abline(h=0)
r2<-result2$residuals
plot(r2,type='o')
abline(h=0)#because of so many facts we get ,the picture is not clear enough to judge wheter the residuals are independent or not.
Box.test(result1$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.3087
Box.test(result2$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.9601
#Judging from the p-value ,ar(1),ar(2)are both accepted.so we should choose ar(1) instead of ar(2) which regards as the fitted model just judging from aic
#5、 考察不同估计方法得到的参数估计结果，即输入如下语句：
arima(x,order=c(1,0,0),method="ML")
#Coefficients:
# ar1 intercept
# -0.6565 0.0026
#s.e. 0.0075 0.0060
#sigma^2 estimated as 0.9901
arima(x,order=c(1,0,0),method="CSS")
#Coefficients:
# ar1 intercept
# -0.6566 0.0026
#s.e. 0.0075 0.0060
#sigma^2 estimated as 0.9902
arima(x,order=c(2,0,0),method="ML")
#Coefficients:
# ar1 ar2 intercept
# -0.6666 -0.0154 0.0026
#s.e. 0.0100 0.0100 0.0059
#sigma^2 estimated as 0.9899
arima(x,order=c(2,0,0),method="CSS")
#Coefficients:
# ar1 ar2 intercept
# -0.6667 -0.0154 0.0026
#s.e. 0.0100 0.0100 0.0059
#sigma^2 estimated as 0.9901
#6、 将估计结果记录下来，可以发现三种方法的估计值很接近，均相当接近于理论真值。（这是大样本，因而估计效果较好。）
arima(x,order=c(1,0,0),method="ML")
#Coefficients:
# ar1 intercept
# -0.6565 0.0026
#s.e. 0.0075 0.0060
#sigma^2 estimated as 0.9901
arima(x,order=c(1,0,0),method="CSS")
#Coefficients:
# ar1 intercept
# -0.6566 0.0026
#s.e. 0.0075 0.0060
#sigma^2 estimated as 0.9902
arima(x,order=c(2,0,0),method="ML")
#Coefficients:
# ar1 ar2 intercept
# -0.6666 -0.0154 0.0026
#s.e. 0.0100 0.0100 0.0059
#sigma^2 estimated as 0.9899
arima(x,order=c(2,0,0),method="CSS")
#Coefficients:
# ar1 ar2 intercept
# -0.6667 -0.0154 0.0026
#s.e. 0.0100 0.0100 0.0059
#sigma^2 estimated as 0.9901
#7、 对序列（2），分别用AR(1)模型、AR(2)模型和AR(3)模型对观测数据进行拟合，利用AIC值确定拟合模型的阶数并与模型实际阶数相比较。
x<-ar2.sim
#平稳性检验
par(mfrow=c(2,1))
acf(x)
pacf(x)
library(tseries)
adf.test(x)# p-value = 0.01， stationary
#数据拟合
result1<-arima(x,order=c(1,0,0))
result2<-arima(x,order=c(2,0,0))
result3<-arima(x,order=c(3,0,0))
result11<-arima(x,order=c(1,0,0),include.mean=F)
result21<-arima(x,order=c(2,0,0),include.mean=F)
result31<-arima(x,order=c(2,0,0),include.mean=F)
#记录上述两个模型的AIC值，确定拟合模型的阶数并与模型实际阶数相比较。
AIC(result1)#29244.69
AIC(result2)#28386.86
AIC(result3)#28388.54
AIC(result11)#29245.83
AIC(result21)#28386.59
AIC(result31)#28386.59
#定阶为ar（2），与实际相符
#对残差进行自相关检验
par(mfrow=c(3,1))
r<-result1$residuals
plot(r,type='o')
abline(h=0)
r<-result2$residuals
plot(r,type='o')
abline(h=0)
r2<-result3$residuals
plot(r2,type='o')
abline(h=0)
Box.test(result1$residuals, lag = 1, type = "Ljung-Box")# p-value < 2.2e-16
Box.test(result2$residuals, lag = 1, type = "Ljung-Box")# p-value = 0.8726
Box.test(result3$residuals, lag = 1, type = "Ljung-Box")# p-value = 0.9954
#ar（2）被接受
# 考察不同估计方法得到的参数估计结果，即输入如下语句：
arima(x,order=c(1,0,0),method="ML")
#Coefficients:
# ar1 intercept
# 0.7129 -0.0645
#s.e. 0.0070 0.0363
#sigma^2 estimated as 1.09
arima(x,order=c(1,0,0),method="CSS")
#Coefficients:
# ar1 intercept
# 0.7129 -0.0645
#s.e. 0.0070 0.0364
#sigma^2 estimated as 1.09
arima(x,order=c(2,0,0),method="ML")
#Coefficients:
# ar1 ar2 intercept
# 0.5082 0.2870 -0.0644
#s.e. 0.0096 0.0096 0.0488
#sigma^2 estimated as 0.9999
arima(x,order=c(2,0,0),method="CSS")
#Coefficients:
# ar1 ar2 intercept
# 0.5082 0.2871 -0.0656
#s.e. 0.0096 0.0096 0.0489
#sigma^2 estimated as 1
arima(x,order=c(3,0,0),method="ML")
#Coefficients:
# ar1 ar2 ar3 intercept
# 0.5066 0.2842 0.0056 -0.0644
#s.e. 0.0100 0.0108 0.0100 0.0491
#sigma^2 estimated as 0.9999
arima(x,order=c(3,0,0),method="CSS")
#Coefficients:
# ar1 ar2 ar3 intercept
# 0.5066 0.2842 0.0056 -0.0654
#s.e. 0.0100 0.0108 0.0100 0.0491
#sigma^2 estimated as 1
#8、 同样对序列（3）、（4）、（5）进行如上操作，分别取拟合模型为MA(1)、MA(2)，MA(1)、MA(2)、MA(3)及ARMA(1,1)、ARMA(2,1)，ARMA(1,2)。将结果记录下来。
#######################对模型3###########################
x<-ma3.sim
#平稳性检验(略）
#数据拟合
result1<-arima(x,order=c(0,0,1))
result2<-arima(x,order=c(0,0,2))
result11<-arima(x,order=c(0,0,1),include.mean=F)
result21<-arima(x,order=c(0,0,2),include.mean=F)
#记录上述两个模型的AIC值，确定拟合模型的阶数并与模型实际阶数相比较。
AIC(result1)#28533.94
AIC(result2)#28534.67
AIC(result11)#28531.99
AIC(result21)#28532.72
#定阶为ma（2），与实际不相符
#对残差进行自相关检验
Box.test(result1$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.5871
Box.test(result2$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.985
#ma（1）被接受,且aic的值相差不大，故选择ma(1)是合理的
# 考察不同估计方法得到的参数估计结果，即输入如下语句：
arima(x,order=c(0,0,1),method="ML")
#Coefficients:
# ma1 intercept
# -0.4851 -0.0012
#s.e. 0.0088 0.0052
#sigma^2 estimated as 1.015
arima(x,order=c(0,0,1),method="CSS")
#Coefficients:
# ma1 intercept
# -0.4851 -0.0012
#s.e. 0.0088 0.0052
#sigma^2 estimated as 1.015
arima(x,order=c(0,0,2),method="ML")
#Coefficients:
# ma1 ma2 intercept
# -0.4793 -0.0113 -0.0012
#s.e. 0.0101 0.0100 0.0051
#sigma^2 estimated as 1.015
arima(x,order=c(0,0,2),method="CSS")
#Coefficients:
# ma1 ma2 intercept
# -0.4794 -0.0113 -0.0012
#s.e. 0.0101 0.0100 0.0051
#sigma^2 estimated as 1.015
#######################对模型4###########################
x<-ma4.sim
#平稳性检验（略）
#数据拟合
result1<-arima(x,order=c(0,0,1))
result2<-arima(x,order=c(0,0,2))
result3<-arima(x,order=c(0,0,3))
result11<-arima(x,order=c(0,0,1),include.mean=F)
result21<-arima(x,order=c(0,0,2),include.mean=F)
result31<-arima(x,order=c(0,0,3),include.mean=F)
#记录上述两个模型的AIC值，确定拟合模型的阶数并与模型实际阶数相比较。
AIC(result1)#28808.24
AIC(result2)#28232.16
AIC(result3)#28233.69
AIC(result11)#28806.27
AIC(result21)#28230.23
AIC(result31)#28231.77
#定阶为ma（2），与实际相符
#对残差进行自相关检验
Box.test(result1$residuals, lag = 1, type = "Ljung-Box")#p-value < 2.2e-16
Box.test(result2$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.8725
Box.test(result3$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.9924
#ma（2）被接受
# 考察不同估计方法得到的参数估计结果，即输入如下语句：
arima(x,order=c(0,0,1),method="ML")
#Coefficients:
# ma1 intercept
# -0.8545 -0.0003
#s.e. 0.0070 0.0015
#sigma^2 estimated as 1.043:
arima(x,order=c(0,0,1),method="CSS")
#Coefficients:
# ma1 intercept
# -0.8547 -0.0003
#s.e. 0.0070 0.0015
#sigma^2 estimated as 1.043
arima(x,order=c(0,0,2),method="ML")
#Coefficients:
# ma1 ma2 intercept
# -0.6473 -0.2424 -0.0003
#s.e. 0.0096 0.0096 0.0011
#sigma^2 estimated as 0.9845
arima(x,order=c(0,0,2),method="CSS")
#Coefficients:
# ma1 ma2 intercept
# -0.6473 -0.2426 -0.0003
#s.e. 0.0096 0.0096 0.0011
#sigma^2 estimated as 0.9846
arima(x,order=c(0,0,3),method="ML")
#Coefficients:
# ma1 ma2 ma3 intercept
# -0.6456 -0.2381 -0.0070 -0.0003
#s.e. 0.0100 0.0116 0.0102 0.0011
#sigma^2 estimated as 0.9845:
arima(x,order=c(0,0,3),method="CSS")
#Coefficients:
# ma1 ma2 ma3 intercept
# -0.6455 -0.2381 -0.0072 -0.0003
#s.e. 0.0100 0.0116 0.0102 0.0011
#sigma^2 estimated as 0.9846
############################模型5##############################
x<-arma5.sim
#平稳性检验（略）
#数据拟合
result1<-arima(x,order=c(1,0,1))
result2<-arima(x,order=c(2,0,1))
result3<-arima(x,order=c(1,0,2))
result11<-arima(x,order=c(1,0,1),include.mean=F)
result21<-arima(x,order=c(2,0,1),include.mean=F)
result31<-arima(x,order=c(1,0,2),include.mean=F)
#记录上述两个模型的AIC值，确定拟合模型的阶数并与模型实际阶数相比较。
AIC(result1)#28384.23
AIC(result2)#28384.11
AIC(result3)#28384.06
AIC(result11)#28383.57
AIC(result21)#28383.4
AIC(result31)#28383.36
#定阶为arma（1，2），与实际相符
#对残差进行自相关检验
Box.test(result1$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.6296
Box.test(result2$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.9927
Box.test(result3$residuals, lag = 1, type = "Ljung-Box")#p-value = 0.9979
#arma（1,2）被接受，其实arma（1,1）也是可以被接受的。
# 考察不同估计方法得到的参数估计结果，即输入如下语句：
arima(x,order=c(1,0,1),method="ML")
#Coefficients:
# ar1 ma1 intercept
# 0.8102 -0.3983 0.0367
#s.e. 0.0098 0.0154 0.0317
#sigma^2 estimated as 0.9997
arima(x,order=c(1,0,1),method="CSS")
#Coefficients:
# ar1 ma1 intercept
# 0.8102 -0.3982 0.0359
#s.e. 0.0098 0.0154 0.0317
#sigma^2 estimated as 0.9997
arima(x,order=c(1,0,2),method="ML")
#Coefficients:
# ar1 ma1 ma2 intercept
# 0.8202 -0.4033 -0.0180 0.0367
#s.e. 0.0115 0.0153 0.0122 0.0322
#sigma^2 estimated as 0.9995
arima(x,order=c(1,0,2),method="CSS")
#Coefficients:
# ar1 ma1 ma2 intercept
# 0.8200 -0.4032 -0.0178 0.0360
#s.e. 0.0115 0.0153 0.0122 0.0322
#sigma^2 estimated as 0.9995
arima(x,order=c(2,0,1),method="ML")
#Coefficients:
# ar1 ar2 ma1 intercept
# 0.8640 -0.0357 -0.4473 0.0367
#s.e. 0.0375 0.0244 0.0358 0.0322
#sigma^2 estimated as 0.9995
arima(x,order=c(2,0,1),method="CSS")
#Coefficients:
# ar1 ar2 ma1 intercept
# 0.8612 -0.0340 -0.4445 0.0358
#s.e. 0.0374 0.0244 0.0358 0.0321
#sigma^2 estimated as 0.9997
#将本次上机实验的结果进行整理，得出你的结论并写出实验报告。
#1、时间序列验证AIC准则定阶的正确性；如果aic的差异大时，用aic定阶比较合理。如果差异不大，那么考虑对残差进行Box.test是比较合理的.我们不能盯住aic不放，结合acf与pacf也是个好办法
#2、利用样本自相关函数、样本偏自相关函数选择模型；
#3、用不同方法估计模型参数。对于大样本，两者的差别不大，但是用阶估计对于含有ma的模型，就十分不好。对于ma模型，相差一阶的差别也不是十分大
扫一扫
3088
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
