Pasted

# ==================================================================================
# r-script para testar ajuste de modelos PAR(p)
# input: Energias medias mensais historicas para REEs (EAF_NE.csv, ...)
# output: relatorio com ordens PAR(p) escolhida via AIC e BIC
# Ref: Estacionaridade das energias afluentes aos REEs
# ==================================================================================
#
# ***********************************************************************************************
# Clean variables & Set working directory
rm(list=ls());                     # limpa vari?veis
# workdir <- "D:/ProjetosONS/Cursos/Alex/R-project/EAFs/2018/";          # diretorio de trabalho
workdir <- file.path("D:", "ProjetosONS", "Cursos", "Alex", "R-project", "EAFs", "2018") 
# workdir <- file.path("/home", "joari", "ProjetosONS", "Cursos", "Alex", "R-project", "ENAs", "2018", "readtxt") 
setwd(workdir);
# ***********************************************************************************************
source("functions.R")
source("readEAFh.R")
#
# ***********************************************************************************************
#
criterion <- "bic"
for (iREE in 1:totREE) {
# for (iREE in 8:8) {
  print(paste0("iREE = ", iREE))
  EAF <- EAFh[iREE,,2:13]
  arquivo <- paste("EAF-", NameREE[[iREE]], "-", toupper(criterion), sep = "")
  saida <- paste(arquivo, ".out", sep = "")
  saida_Full <- paste(arquivo, "-Full.out", sep = "")
  #
  # ***********************************************************************************************
  #
  # fit parp (pear package), using library("pear")
  maxEAF <- 1.05*max(EAF)
  savePARp <- list()
  cat(paste0("\n ", toupper(criterion), " model order\n", sep=""), file=saida, append=FALSE)
  cat(paste0("\n ", toupper(criterion), " model order parameteres\n", sep=""), file=saida_Full, append=FALSE)
  title <- paste(" Anoini Anofin   Total", paste(formatC(month.abb,format="s",width=6), collapse = " "), sep = "")
  cat(title, file=saida, append=TRUE)
  cases <- seq(1, 66, by = 5)
  # for (inicio in cases){
  for (i in 1:length(cases)){
    anoini <- anoIni + cases[i] - 1
    out <- paste("\n")
    out <- paste(out, 
                 formatC(anoini, digits=0, format="f", width=6),
                 formatC(anoFim, digits=0, format="f", width=6),
                 formatC(anoFim - anoini + 1, digits=0, format="f", width=7), sep = " ")
    cat(out, file=saida, append=TRUE)
    cat(out, file=saida_Full, append=TRUE)
    tmp <- EAF[cases[i]:totYears,]
    plot_EAFh (tmp, maxEAF, anoIni, anoFim, arquivo, anoini)
    # padronizando tmp
    EAFpad <- scale(tmp)
    savePARp[[i]] <- try(ajusta_parp (EAFpad, anoFim, criterion, saida, saida_Full), silent = TRUE)
  }
  #
  aa   <- array(0, dim = c(length(cases), 12))
  ymed <- array(0, dim = c(length(cases), 12))
  ysup <- array(0, dim = c(length(cases), 12))
  yinf <- array(0, dim = c(length(cases), 12))
  minY <- 1.0
  for (mes in 1:12) {
    # get identified model order for each month (mes) for every sample size
    for (i in seq(1,length(cases))) {
      try({
        if (is.numeric(savePARp[[i]][[1]][[mes]])) {aa[i, mes] <- savePARp[[i]][[1]][[mes]]}
      })
    }
    cat(paste0("mes ", formatC(mes, digits=NULL,format="d",width=3)))
    cat(formatC(aa[, mes], digits=NULL,format="d",width=3), sep="")
    cat("\n")
    if (all(aa[, mes] == 1, na.rm = TRUE)) {
      print(paste0(month.abb[mes], " ordem 1 para todos os casos"))
      for (i in seq(1,length(cases))) {
        ymed[i, mes] <- savePARp[[i]][[2]][mes,1]
        delta <- 1.96*savePARp[[i]][[3]][mes,1]/sqrt(anoFim - anoIni - cases[i] + 2)
        ysup[i, mes] <- ymed[i, mes] + delta
        yinf[i, mes] <- ymed[i, mes] - delta
        minY <- min(minY, yinf[i, mes])
      }
    }
  }
  # Barplot for all modes with order == 1 for every month
  anoI = (anoIni + cases - 1)
  maxY <- 1
  for (mes in 1:12) {
    if (all(aa[, mes] == 1)) {
      png(file=paste(arquivo, " - ", month.abb[mes],".png",sep=""));
      phi_barplot <- barplot(ymed[, mes], names.arg = anoI, col=c("blue" , "skyblue") ,
                             ylim = c(minY, maxY), xlab = "ano inicial", ylab = "Model order",
                             main = paste(arquivo, " - ", month.abb[mes], sep = ""),
                             xpd=FALSE)
      axis(side = 1,at = phi_barplot, labels = FALSE)
      arrows(phi_barplot, ysup[, mes], phi_barplot, yinf[, mes], angle = 90, code = 3, length = 0.1)
      abline(h = 1, col = "blue", lwd = 2)
      box(bty = "l")
      dev.off();
    }
  }
  plot_order(aa, cases, anoI, arquivo)
}