Como está começando, segue aqui uma sugestao valiosa: leia o "posting guide" da lista em resumo, como diz um mui solicito participante da lista: lembre de "ajudar a quem quer te ajudar" "Tenho tido problemas" é muito vago e seu codigo muito longo, provavelmente com trechos desnecessários para a postagem. Entenda o que é um CMR (código mínimo repreduzível). Seja claro, minimo e lembre sempre que quem se dispoe a ler e ajudar está usando tempo para isto, portanto use o seu também para fazer uma postagem objetiva Boas postagens! P.J. On Fri, 30 Mar 2012, Vinicius Macedo Magalhães wrote:
Fala, galera. Primeira vez usando este fórum. Tenho tido problemas para rodar meu código para previsão sazonal via ARIMA. Abaixo segue parte do código.
s.ts Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 2004 118 143 169 158 143 135 135 140 135 125 120 120 2005 143 158 180 180 150 150 153 148 150 145 145 143 2006 165 180 190 168 163 163 168 175 180 160 155 150 2007 152 163 185 163 158 175 175 178 168 168 160 160 2008 155 170 195 195 195 185 185 165 165 153 154 159 ##### Início de ARIMA ##### linhaX <- 3 dados <- matrix() nome <- names(db) i=4 get.best.arima <- function (x.ts, maxord=c(1,1,1,1,1,1)) { best.aic <- 1e8 n <- length(x.ts) for (p in 0:maxord[1]) for (d in 0:maxord[2]) for (q in 0:maxord[3]) for (P in 0:maxord[4]) for (D in 0:maxord[5]) for (Q in 0:maxord[6]) { fit <- arima(x.ts, order=c(p,d,q),seasonal=list(order=c(P,D,Q),period=12), method = "CSS") fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef) if (fit.aic < best.aic) { best.aic <- fit.aic best.fit <- fit best.model <- c(p,d,q,P,D,Q) } } list(best.aic,best.fit,best.model) }
s.ts <- ts(db[,i],start=2004, fre=12) best.arima.ss <- get.best.arima(log(s.ts),maxord = c(2,2,2,2,2,2)) best.fit.ss <- best.arima.ss[[2]] best.arima.ss[[3]] pdf(paste("AR.", i-4, ".pdf", sep=""),width=7,height=5) ts.plot( cbind( window(s.ts, start=2004),exp(predict(best.fit.ss,24)$pred)),col=c("black","green")) #par(new=TRUE) #ts.plot(s.ts,best.fit.ss,col=c("black","red")) legend("bottomright",c("Real","Previsão_Passado","Previsão_Futuro"),col=1:3,lty=1,cex=0.7,merge = TRUE,bg = 'gray90') dev.off() report<-function(Actual,Forecast,Net,R2) { result <- new.env() if (R2==0) R2<-1-(sum((Forecast-Actual)^2)/sum((Actual-mean(Actual))^2)) result$Rsquare<-R2 result$MAPE <- 100*(sum(abs((Actual-Forecast))/Actual))/nrow(Forecast) ## Mean Absolute Percent Error result$CFE <- sum((Actual-Forecast)) ## Cumulative Forecast Error result$ME <- result$CFE/nrow(Forecast) ## Mean Error result$MSE <- sum((Actual-Forecast)^2)/nrow(Forecast) ## Mean Squared Error result$RMSE <- sqrt(result$MSE) ## Root Mean Squared Error result$MAD <- (sum(abs(Actual-Forecast))/nrow(Forecast)) ## Mean Absolute Deviation result$SUMMARY <- summary(Net) return(result) } ## Fim da função report pred.past<- db[,colY] previsao <- predict(best.fit.ss) forecast <- as.matrix(pred.past[1:length(pred.past)]) reportAR <- report(db[,colY],forecast,best.fit.ss,0) dados <- rbind(dados,"Produto:"=nome[colY],R2=reportAR$Rsquare,MAPE=reportAR$MAPE, CFE=reportAR$CFE,ME=reportAR$ME,MSE=reportAR$MSE,RMSE=reportAR$RMSE,MAD=reportAR$MAD) write.table(dados,sep=",",file="ReportAR.csv",)
} ## Fim do for return(dados) ##### Fim de ARIMA #####
-- Att.,
Vinicius Macedo Magalhães (21) 9584-1533