functions.R

# Functions for this study
# ==================================================================================
# Load libraries
library(pear)
# ==================================================================================
plot_order <- function(aa, cases, anoI, arquivo){
  maxY <- max(aa)
  validSample <- NULL
  for (j in 1:length(cases)) {
    if (aa[j,1] != 0) {
      validSample = c(validSample, (anoIni + cases[j] - 1))
    }
  }
  png(file=paste(arquivo, " - Model order",".png", sep=""));
  par(mar = par()$mar + c(0,0,0,4))
  if (aa[1,1] != 0){
    plot(1:12, y = aa[1,], type = "o", pch=1,
         xlab = "Month", ylab = "Model order",
         ylim = c(0, maxY), main = arquivo, xaxt='n')
    axis(side = 1, at = 1:12, labels=month.abb)
  }
  for (i in 2:length(cases)) {
    if (aa[i,1] != 0) {
      lines(aa[i,], type = "o", pch=i, col=rainbow(14)[i])
    }
  }
  # ll <- seq(1:14)
  # legend("topleft", inset=.05, legend = anoI, title="Sample", pch = ll, bty = "n", fill=rainbow(14))#, horiz=F)
  legend("right", inset = c(-0.21,0), legend = validSample, xpd = TRUE,
         horiz = FALSE, col = rainbow(14), lty = 1, bty = "n")
  par(mar=c(5, 4, 4, 2) + 0.1)
  dev.off()
}
# ==================================================================================
plot_EAFh <-  function(EAF, maxEAF, anoIni, anoFim, arquivo, anoini){
  tmp <- as.vector(t(EAF))
  EAF_ts <- ts(tmp, end = anoFim, frequency = 12)

  png(file=paste(arquivo, "-TS-inic-", anoini, ".png", sep=""));
  plot(EAF_ts, xlim = c(anoIni,anoFim), ylim = c(0, maxEAF))
  title(main = arquivo)
  mtext(paste0("Ano inicial = ", anoini), side = 3, line = 0.5)
  #summary(EAF_ts)
  #hist(EAF_ts,prob=TRUE)
  dev.off()

  png(file=paste(arquivo, "-BP-inic-", anoini, ".png", sep=""));
  boxplot(EAF_ts ~ cycle(EAF_ts), names=month.abb, ylim = c(0, maxEAF))#, main = arquivo)
  title(main = arquivo)
  mtext(paste0("Ano inicial = ", anoini), side = 3, line = 0.5)
  dev.off()

  png(file=paste(arquivo, "-PePACF-inic-", anoini, ".png", sep=""));
  outpepacf <- pepacf(EAF_ts, 15, plot=TRUE)
  title(main = arquivo)
  mtext(paste0("Ano inicial = ", anoini), side = 3, line = 0.5)
  dev.off()
}
# ==================================================================================
ajusta_parp <- function(EAF, anoFim, criterion, saida, saida_Full){
  tmp <- as.vector(t(EAF))
  EAF_ts <- ts(tmp, end = anoFim, frequency = 12)
  parpIC <- list()
  #  for (criterion in c("aic", "bic")){
  # for (criterion in c("aic")){
    parp<-pear(EAF_ts,ic=criterion)
    # parp <- pear(EAF_ts,ic="aic")
    parpIC <- c(parpIC , parp)
    #parp$model.orders
    # IC <- toupper(criterion)
    cat(paste(formatC(parp$model.orders,digits=NULL,format="d",width=6), sep=""), file=saida, append=TRUE)
    title <- paste("\n", toupper(criterion)," phi coefficients\n", sep = "" )
    prt.phi(title, parp$phi, parp$model.orders, saida_Full)
    title <- paste(toupper(criterion)," phi standard deviation \n", sep = "" )
    prt.phi(title, parp$se.phi, parp$model.orders, saida_Full)
  # }
  return(parpIC)
}
# ==================================================================================
# print PAR(p) phi or sd.phi coefficients as formatted table
prt.phi <- function(title, phi, model.orders, saida){
  cat(title,file=saida,append=T);
  ordMax <- max(model.orders)
  cat(" month", paste(formatC("lag_",format="s",width=11), formatC(1:ordMax,format="d",flag="0",width=2),sep=""),
      file=saida, fill=13*22+6, append=T);
  for (i in 1:12){
    cat(formatC(phi[i,1:model.orders[[i]]],digits=9,format="f",width=13), labels=paste(formatC(month.abb[i],format="s",width=6)),
        file=saida, fill=13*22+6, append=T);
  }
}
# ==================================================================================
#