2

Preciso verificar o percentual de vezes que AIC e BIC escolhem o verdadeiro modelo. Para tanto, terá que realizar um experimento de Monte Carlo.
Especificamente, devem ser gerados 1000 processos AR(2) e ARMA(1,1). Os resultados de cada modelo devem ser comparados com os de modelos alternativos.
Para facilitar, utilize os seguintes modelos alternativos para comparar comparar com o AR(2): ARMA(2,1), ARMA(1,1), AR(1,0), ARMA(1,2). Para comparar com o ARMA(1,1), considere as seguintes possibilidades: ARMA(1,0), ARMA(2,0), ARMA(0,2) e ARMA(1,2). O número de observações do processo yt é igual a 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))

} 

Está dando o seguinte erro:

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

Como corrigir este erro?

3
  • Marta, claramente isto é um exercício. Simplesmente apontar pra você onde está o erro não te ajudará a aprender sobre R. Procure rodar o código passo a passo e observar o resultado de cada etapa. Recomendo também que você se atente à formatação do código, ao use o operador <- para definição de objetos ao invés do símbolo de igual e procure utilizar as funções da família *apply ao invés de tantos for loops.
    – Molx
    13/04/2016 às 4:07
  • Não vejo muito sentido nesse posicionamento, @Molx, o erro aqui é um erro bastante comum e certamente há aprendizado em identificar esse erro. 14/04/2016 às 16:55
  • 1
    @CarlosCinelli Ao meu ver, quando a pessoa está tendo um erro simples em um código tão extenso é porque está com falhas em fundamentos básicos da linguagem ou muitas vezes nem é o autor do código. Nesses casos simplesmente dar a resposta não resolve nada, porque este tipo de erro é simples de ser resolvido se a pessoa aprende outras coisas antes. Me parece que etapas do aprendizado foram puladas e nada será aprendido, apenas a resposta copiada, e não será útil para mais ninguém. Basta ver que o erro não tem nada a ver com o título ou a descrição da pergunta. Por isso, minhas sugestões.
    – Molx
    14/04/2016 às 20:56

1 Resposta 1

1

O erro na verdade é simples, é um erro de indexação no loop.

Veja que no primeiro loop, nas linhas em que você atualiza as matrizes, você está usando rep ao invés de t.

Isto é, nas linhas:

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)

Você está sempre atualizando a última linha pois você está passando a variável rep. Basta trocar rep por t que as matrizes serão preenchidas corretamente.

Segue o primeiro loop corrigido:

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

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .