2
I need to check the percentage of times AIC and BIC choose
the real model. For this, you will have to carry out a Monte Carlo experiment.
Specifically, 1000 AR(2) and ARMA(1,1) processes must be generated. The
results of each model shall be compared with those of alternative models.
For ease, use the following alternative models to compare with
o AR(2): ARMA(2,1), ARMA(1,1), AR(1,0), ARMA(1,2). To compare with
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 yt process
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’re making the following mistake:
Error in min.aic.ar2[t] = which(aic.ar2[t, ] == min(aic.ar2[t, ])) :
replacement has zero length
How to correct this error?
Marta, clearly this is a drill. Simply pointing to you where the error is will not help you learn about R. Try to run the code step by step and observe the outcome of each step. I also recommend that you stick to code formatting, when using the operator
<-for setting objects instead of the equals symbol and try to use the family functions*applyinstead of so many loops.– Molx
I don’t see much point in this positioning, @Molx, the error here is a very common mistake and certainly there is learning to identify this error.
– Carlos Cinelli
@Carloscinelli The way I see it, when a person is having a simple error in such a large code is because they are lacking in basic language fundamentals or are often not even the author of the code. In such cases simply giving the answer does not solve anything, because this kind of mistake is simple to be solved if one learns other things first. It seems to me that stages of learning have been skipped and nothing will be learned, only the response copied, and will be of no use to anyone else. Just note that the error has nothing to do with the title or description of the question. So my suggestions.
– Molx