[R-br] Ajuda - Arima Sazonal

[Vinicius Macedo Magalhães]

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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120330/92eae79d/attachment.html>