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?