information criterion Akaike time series

2

I need to check the percentage of times AIC and BIC choose the true model. To do so, you will have to perform a Monte Carlo experiment. Specifically, 1000 AR (2) and ARMA (1,1) processes must be generated. The results of each model should be compared with those of alternative models. To facilitate, use the following alternative models to compare the AR (2): ARMA (2,1), ARMA (1,1), AR (1,0), ARMA (1,2). To compare with the ARMA (1,1), consider the following possibilities: ARMA (1,0), ARMA (2,0), ARMA (0.2) and ARMA (1,2). The number of observations of the process yt  is equal to 500.

rep=200
nmodels=5
#Inicializa matrizes para armazenar informações.
aic.ar2=matrix(,nrow=rep,ncol=nmodels)
bic.ar2=matrix(,nrow=rep,ncol=nmodels)
aic.arma11=matrix(,nrow=rep,ncol=nmodels)
bic.arma11=matrix(,nrow=rep,ncol=nmodels)
colnames(aic.ar2)<-c("ar2","arma(2,1)","arma(1,1)","ar1","arma(1,2)")
colnames(aic.arma11)<-c("arma(1,1)","arma(1,0)","arma(2,0)","arma(0,2)","arma(1,2)")

#II. Loop
#i) Gera uma sequência de AR(2) e uma de ARMA(1,1).
 for (t in 1: rep){
 y.ar=arima.sim(list(ar=c(0.4,-0.3)),n=190)
 y.arma=arima.sim(list(order=c(1,0,1),ar=0.7, ma=.5),n=190)

#ii) Estimar AR(2) e modelos para comparar com AR(2)
 mod1=arima(y.ar,order=c(2,0,0))
 mod2=arima(y.ar,order=c(2,0,1))
 mod3=arima(y.ar,order=c(1,0,1))
 mod4=arima(y.ar,order=c(1,0,0))
 mod5=arima(y.ar,order=c(1,0,2))

#i2) Estimar ARMA(1,1) e modelos para comparar com ARMA(1,1)
 mod1a=arima(y.arma,order=c(1,0,1))
 mod2a=arima(y.arma,order=c(1,0,0))
 mod3a=arima(y.arma,order=c(2,0,0))
 mod4a=arima(y.arma,order=c(0,0,2))
 mod5a=arima(y.arma,order=c(1,0,2))

#iii) Guardar os valores do AIC
 aic=c(mod1$aic,mod2$aic,mod3$aic,mod4$aic,mod5$aic)
 aic.ar2[rep, 1:nmodels] = aic

 bic1=AIC(mod1,k=log(length(y.ar)))
 bic2=AIC(mod2,k=log(length(y.ar)))
 bic3=AIC(mod3,k=log(length(y.ar)))
 bic4=AIC(mod4,k=log(length(y.ar)))
 bic5=AIC(mod5,k=log(length(y.ar)))
 bic.ar2[rep,1:nmodels]=c(bic1,bic2,bic3,bic4,bic5)

 aic=c(mod1a$aic,mod2a$aic,mod3a$aic,mod4a$aic,mod5a$aic)
 aic.arma11[rep,1:nmodels] = aic

 bic1a=AIC(mod1a,k=log(length(y.arma)))
 bic2a=AIC(mod2a,k=log(length(y.arma)))
 bic3a=AIC(mod3a,k=log(length(y.arma)))
 bic4a=AIC(mod4a,k=log(length(y.arma)))
 bic5a=AIC(mod5a,k=log(length(y.arma)))
 bic.arma11[rep,1:nmodels]=c(bic1a,bic2a,bic3a,bic4a,bic5a)

 } #fecha loop

#III. Comparação
min.aic.ar2=1:rep
min.bic.ar2=1:rep
min.aic.arma11=1:rep
min.bic.arma11=1:rep

for (t in 1:rep){
 min.aic.ar2[t]=which(aic.ar2[t,]==min(aic.ar2[t,]))
 min.aic.arma11[t]=which(aic.arma11[t,]==min(aic.arma11[t,]))

 min.bic.ar2[t]=which(bic.ar2[t,]==min(bic.ar2[t,]))
 min.bic.arma11[t]=which(bic.arma11[t,]==min(bic.arma11[t,]))

 } 

best.aic.ar2=matrix(0,ncol=5,nrow=1)
best.bic.ar2=matrix(0,ncol=5,nrow=1)
colnames(best.aic.ar2)<-c("ar2"," arma(2,1)"," arma(1,1)"," ar1"," arma(1,2)")
colnames(best.bic.ar2)<-c("ar2"," arma(2,1)"," arma(1,1)"," ar1"," arma(1,2)")

best.aic.arma11=matrix(0,ncol=5,nrow=1)
best.bic.arma11=matrix(0,ncol=5,nrow=1)
colnames(best.aic.arma11)<-c("arma(1,1)","arma(1,0)","arma(2,0)","arma(0,2)","arma(1,2)")
colnames(best.bic.arma11)<-c("arma(1,1)","arma(1,0)","arma(2,0)","arma(0,2)","arma(1,2)")

for (i in 1:nmodels){
 best.aic.ar2[,i]=length(which(min.aic.ar2==i))
 best.bic.ar2[,i]=length(which(min.bic.ar2==i))

 best.aic.arma11[,i]=length(which(min.aic.arma11==i))
 best.bic.arma11[,i]=length(which(min.bic.arma11==i))

} 

You are giving the following error:

Error in min.aic.ar2[t] = which(aic.ar2[t, ] == min(aic.ar2[t, ])) : 
  replacement has zero length

How to fix this error?

    
asked by anonymous 13.04.2016 / 01:46

1 answer

1

The error is actually simple, it is an indexing error in the loop.

Notice that in the first loop, in the rows where you update the arrays, you are using rep instead of t .

That is, in the lines:

aic.ar2[rep, 1:nmodels] = aic

bic.ar2[rep,1:nmodels]=c(bic1,bic2,bic3,bic4,bic5)

aic.arma11[rep,1:nmodels] = aic

bic.arma11[rep,1:nmodels]=c(bic1a,bic2a,bic3a,bic4a,bic5a)

You are always updating the last line because you are passing the variable rep . It is enough to change rep by t that the matrices will be filled correctly.

The first corrected loop follows:

rep=200
nmodels=5

#Inicializa matrizes para armazenar informações.
aic.ar2=matrix(,nrow=rep,ncol=nmodels)
bic.ar2=matrix(,nrow=rep,ncol=nmodels)
aic.arma11=matrix(,nrow=rep,ncol=nmodels)
bic.arma11=matrix(,nrow=rep,ncol=nmodels)
colnames(aic.ar2)<-c("ar2","arma(2,1)","arma(1,1)","ar1","arma(1,2)")
colnames(aic.arma11)<-c("arma(1,1)","arma(1,0)","arma(2,0)","arma(0,2)","arma(1,2)")
aic.ar2
t = 1
#II. Loop
#i) Gera uma sequência de AR(2) e uma de ARMA(1,1).
for (t in 1: rep){
  y.ar=arima.sim(list(ar=c(0.4,-0.3)),n=190)
  y.arma=arima.sim(list(order=c(1,0,1),ar=0.7, ma=.5),n=190)

  #ii) Estimar AR(2) e modelos para comparar com AR(2)
  mod1=arima(y.ar,order=c(2,0,0))
  mod2=arima(y.ar,order=c(2,0,1))
  mod3=arima(y.ar,order=c(1,0,1))
  mod4=arima(y.ar,order=c(1,0,0))
  mod5=arima(y.ar,order=c(1,0,2))

  #i2) Estimar ARMA(1,1) e modelos para comparar com ARMA(1,1)
  mod1a=arima(y.arma,order=c(1,0,1))
  mod2a=arima(y.arma,order=c(1,0,0))
  mod3a=arima(y.arma,order=c(2,0,0))
  mod4a=arima(y.arma,order=c(0,0,2))
  mod5a=arima(y.arma,order=c(1,0,2))

  #iii) Guardar os valores do AIC
  aic=c(mod1$aic,mod2$aic,mod3$aic,mod4$aic,mod5$aic)
  aic.ar2[t, 1:nmodels] = aic

  bic1=AIC(mod1,k=log(length(y.ar)))
  bic2=AIC(mod2,k=log(length(y.ar)))
  bic3=AIC(mod3,k=log(length(y.ar)))
  bic4=AIC(mod4,k=log(length(y.ar)))
  bic5=AIC(mod5,k=log(length(y.ar)))
  bic.ar2[t,1:nmodels]=c(bic1,bic2,bic3,bic4,bic5)

  aic=c(mod1a$aic,mod2a$aic,mod3a$aic,mod4a$aic,mod5a$aic)
  aic.arma11[t,1:nmodels] = aic

  bic1a=AIC(mod1a,k=log(length(y.arma)))
  bic2a=AIC(mod2a,k=log(length(y.arma)))
  bic3a=AIC(mod3a,k=log(length(y.arma)))
  bic4a=AIC(mod4a,k=log(length(y.arma)))
  bic5a=AIC(mod5a,k=log(length(y.arma)))
  bic.arma11[t,1:nmodels]=c(bic1a,bic2a,bic3a,bic4a,bic5a)

} #fecha loop
    
14.04.2016 / 18:53