时间序列作业

 

第二次上机实验：模拟MA模型和ARMA模型

 作业内容请访问：http://www.doc88.com/p-209946619868.html

参考解答：

#2--5上机结果
n<-100;     
x1<-0;
a2<-0;
a3<-0;
x<-0;      
for(i in -50:n)          
{a1<-rnorm(1);          
x1<-a1+0.65*a2+0.24*a3;
a3<-a2;
a2<-a1;                    
if(i>0)
x[i]<-x1}        
x
mean(x)#-0.007158113
var(x)#1.532393
fivenum(x)#-2.33380162 -0.96968524  0.01153527  0.78820663  2.37949770
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture1

par(mfrow=c(1,2))
acf(x)
pacf(x)#see picture 2

#6
set.seed(100)
n<-100;     
x1<-0;
a2<-0;
a3<-0;
x<-0;      
for(i in -50:n)          
{a1<-rnorm(1);          
x1<-a1-0.65*a2-0.24*a3;
a3<-a2;
a2<-a1;                    
if(i>0)
x[i]<-x1}        
x
mean(x)#0.009997361
var(x)# 1.745586
fivenum(x)#-3.54820575 -0.91206753 -0.02810196  0.85317906  3.06819255
par(mfrow=c(1,1))
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture3

par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture 4

#7
set.seed(100)
n<-100;     
x1<-0;
a2<-0;
a3<-0;
x<-0;      
for(i in -50:n)          
{a1<-rnorm(1);          
x1<-a1-0.65*a2+0.24*a3;
a3<-a2;
a2<-a1;                    
if(i>0)
x[i]<-x1}        
x;
mean(x)#-0.04177703
var(x)# 1.864563
fivenum(x)#-3.46235063 -1.03775369 -0.06152364  0.93842732  2.96942047
par(mfrow=c(1,1))
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture 5

par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture 6

#8
set.seed(100)
x<-arima.sim(n = 100, list( ma = c(0.65, -0.24)))
mean(x)#-0.03727247
var(x)# 1.401805
fivenum(x)#-3.54417536 -0.69506541 -0.04359503  0.66448787  3.33484756
par(mfrow=c(1,1))
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture 7

par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture 8

#8 use function in R directly
set.seed(100)
x<-arima.sim(n = 100, list(ar=c( ma = c(0.65, -0.24)))
mean(x)#-0.03727247
var(x)# 1.401805
fivenum(x)#-3.54417536 -0.69506541 -0.04359503  0.66448787  3.33484756
par(mfrow=c(1,1))
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture 7

par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture 8

#9
set.seed(100)
x<-arima.sim(n = 100, list(ar=c(-0.75,-0.5), ma = c(0.3, 0.4)))+5
mean(x)#4.98
var(x)# 1.692439
fivenum(x)# 1.658756 4.142209 4.996332 5.941509 8.064161
par(mfrow=c(1,1))
plot(x,type="b",xlab="i",ylab=expression(X[i]),col=2)#see picture 9

par(mfrow=c(2,1))
acf(x)
pacf(x)#see picture 10
#10回到graph窗口观察各种序列图形的异同。
#acf均是2步截尾的（ARMA不是）pacf均是拖尾的，arma的acf也是拖尾的
#acf衰减情况不同，有指数衰减的，也有震荡衰减的

#11 模拟其他随机序列，体会不同类型随机序列的样本自相关函数和样本偏自相关函数的差异，尤其是拖尾性与截尾性的区分。
require(graphics)

x1<-arima.sim(n = 63, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
          sd = sqrt(0.1796))
# mildly long-tailed
x2<-arima.sim(n = 63, list(ar=c(0.8897, -0.4858), ma=c(-0.2279, 0.2488)),
          rand.gen = function(n, ...) sqrt(0.1796) * rt(n, df = 5))
par(mfrow=c(2,1))
acf(x1)
pacf(x1)
par(mfrow=c(2,1))
acf(x2)
pacf(x2)
#截尾性与拖尾性其实也没有那么好判断，这与acf的衰减幅度有关。如果衰减幅度特别大的时候，用ARMA拟合效果也不错
#思考：如何利用样本自相关函数和样本偏自相关函数识别序列类型？
#看截尾性与拖尾性，也可以使用广义自相关函数eacf去判断