Pasted

# =============================================================================
# r-script para testar ajuste de modelos PAR(p)
# input: Energias medias mensais historicas para REEs (ENA_NE.csv, ...)
# output: relatorio com ordens PAR(p) escolhida via AIC e BIC
# Ref: Estacionaridade das energias afluentes aos REEs
# =============================================================================
# Remove all objects from the current workspace:
rm(list=ls());
#
# =============================================================================
# Load libraries
library(pear)
#
# =============================================================================
# Functions for this study
ajusta_parp <- function(ENA, anoFim, saida){
  tmp <- as.vector(t(ENA))
  ENA_ts <- ts(tmp, end = anoFim, frequency = 12)
  plot(ENA_ts, xlim = c(1931,2018))
  #summary(ENA_ts)
  #hist(ENA_ts,prob=TRUE)
  #boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
  
  #outpepacf <- pepacf(ENA_ts, 15, plot=TRUE)
  parpIC <- list()
#  for (criterion in c("aic", "bic")){
  for (criterion in c("bic")){
    # parp<-pear(ENA_ts,ic=criterion)
    parp<-pear(ENA_ts,ic="bic")
    parpIC <- c(parpIC , parp)
    #parp$model.orders
    IC <- toupper(criterion)
    #title <- paste("\n", IC," model order\n", sep = "" )
    #cat(title, file=saida, append=TRUE)
    #cat(formatC(month.abb,format="s",width=6), file=saida, append=TRUE)
    #cat(paste("\n"), file=saida, append=TRUE)
    cat(paste(formatC(parp$model.orders,digits=NULL,format="d",width=6), sep=""), file=saida, append=TRUE)
    # title <- paste("\n", IC," phi coefficients\n", sep = "" )
    # prt.phi(title, parp$phi, parp$model.orders, saida)  
    # title <- paste(IC," phi standard deviation \n", sep = "" )
    # prt.phi(title, parp$se.phi, parp$model.orders, saida)
    
  }
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);
  }
}
# =============================================================================
#
# ##############################################################################
# Clean workspace & Set working directory
rm(list=ls());                        # Remove all objects from the current workspace:
# workdir <- "D:/ProjetosONS/Cursos/Alex/R-project/ENAs/2018/"; # diretorio de trabalho
workdir <- "~/ProjetosONS/Cursos/Alex/R-project/ENAs/2018/"; # diretorio de trabalho
setwd(workdir);
#
# ==================================================================================
# le dados de ENAs historicas medias mensais 
arquivo <- "ENA_NE_2016"
# arqin <- paste(arquivo, ".csv", sep = "")
saida <- paste(arquivo, ".out", sep = "")
 
# dado <- read.csv(arqin, header = T, sep = ";") # tabcsv
#
# ==================================================================================
 
# fit parp (pear package) library("pear")
# ENA <- as.matrix(dado[,2:13])
ENA <- Btt[7,,2:13]
objtmp <- list()
# totAnos <- length(ENA[,1]); # comprimento total da amostra
# anoFim <- 2018
cat("\n BIC model order\n", file=saida, append=TRUE)
title <- paste(" Anoini Anofin   Total", paste(formatC(month.abb,format="s",width=6), collapse = " "), sep = "")
cat(title, file=saida, append=TRUE)
for (inicio in seq(1, 51, by = 5)){
# for (inicio in seq(1, 41, by = 40)){
  count <- count + 1
  anoini <- anoIni + inicio - 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)
  tmp <- ENA[inicio:totYears,]
  objtmp <- c(objtmp , ajusta_parp (tmp, anoFim, saida))
}
 
aa <- array()
for (i in seq(1,11)) {
  nObj <- (i-1)*9  + 1
  print(objtmp[[nObj]][[2]])
  aa[i] <- objtmp[[nObj]][[2]]
}
 
print(aa)
# guradar valores do total de anos para x, ou inicio da serie considerada
plot(x = (anoIni + seq(1, 51, by = 5) - 1), y=aa, main=month.abb[2], ylab = "order")
barplot(aa, names.arg = (anoIni + seq(1, 51, by = 5) - 1), 
        xlab = "ano inicial", ylab = "model order", main = paste(arquivo, month.abb[2]))
 
 
seq(1, 41, by = 40)
month.abb[2]
ENA <- as.matrix(dado[,2:13])
totAnos <- length(ENA[,1])
 
tmp <- as.vector(t(ENA))
#tmp <- as.vector(t(ENA[20:85, 1:12]))
 
ENA_ts <- ts(tmp, end = 2018, frequency = 12)
 
 
parp<-pear(ENA_ts,ic="aic")
model <-data.frame()
model <- parp
mydata <- data.frame(parp)
model[1] <- parp
parp$model.orders
parp$phi
parp$se.phi
 
out <- pepacf(ENA_ts, 15, plot=TRUE)
pepacf(ENA_ts, 15, plot=TRUE)
plot(ENA_ts)
summary(ENA_ts)
hist(ENA_ts,prob=TRUE)
boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
acf(as.numeric(ENA_ts))
pacf(as.numeric(ENA_ts))
 
 
 
 
inicio <- 1
anoini <- 1931 + inicio - 1
title <- paste("\n Ano inicial: ", anoini, " Ano final: ", anoFim, " Total : ", anoFim - anoini, sep = "" )
cat(title, file=saida, append=TRUE)
tmp <- ENA[inicio:totAnos,]
ajusta_parp (tmp, anoFim, saida)
  
inicio <- 20
anoini <- 1931 + inicio - 1
title <- paste("\n Ano inicial: ", anoini, " Ano final: ", anoFim, " Total : ", anoFim - anoini, sep = "" )
cat(title, file=saida, append=TRUE)
ENA <- as.matrix(dado[20:totAnos,2:13])
ajusta_parp (ENA, 2018, saida)
 
totAnos <- length(ENA[,1]); # comprimento total da amostra
 
anoFim <- 2016
 
# ==================================================================================
#
# fit parp (pear package) library("pear")
# library(pear)
# ENA_ts = ENA como objeto ts a partir da matriz de dados
# Yts <- ENA[1,]
# for (i in 2:totAnos){
#   Yts <- c(Yts, ENA[i,])
# }
 
tmp <- as.vector(t(ENA))
 
#ENA_ts <- ts(Yts, start = 1931, frequency = 12)
ENA_ts <- ts(tmp, start = 1931, frequency = 12)
 
plot(ENA_ts)
summary(ENA_ts)
hist(ENA_ts,prob=TRUE)
boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
acf(as.numeric(ENA_ts))
pacf(as.numeric(ENA_ts))
 
parp<-pear(ENA_ts,ic="aic")
parp$model.orders
parp$phi
parp$se.phi
 
parp<-pear(ENA_ts, ic="bic")
parp$model.orders
 
# ==================================================================================
 
arquivo <- "ENA_NE_PMO2017"
arqin <- paste(arquivo, ".csv", sep = "")
dado <- read.csv(arqin, header = T, sep = ";") # tabcsv
#
# ==================================================================================
 
ENA <- as.matrix(dado[,2:13])
totAnos <- length(ENA[,1])
 
tmp <- as.vector(t(ENA))
tmp <- as.vector(t(ENA[20:85, 1:12]))
 
ENA_ts <- ts(tmp, end = 2015, frequency = 12)
 
out <- pepacf(ENA_ts, 15, plot=TRUE)
pepacf(ENA_ts, 15, plot=TRUE)
plot(ENA_ts)
summary(ENA_ts)
hist(ENA_ts,prob=TRUE)
boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
acf(as.numeric(ENA_ts))
pacf(as.numeric(ENA_ts))
 
parp<-pear(ENA_ts,ic="aic")
parp$model.orders
parp$phi
parp$se.phi
 
parp<-pear(ENA_ts, ic="bic")
parp$model.orders
 
 
 
 
 
all.equal(parp$model.phi, parp1$model.phi)
 
all(Yts==tmp)
all.equal(Yts,tmp, check.attributes=FALSE)
 
Yts <- ENA[40,]
for (i in 41:totAnos){
  Yts <- c(Yts, ENA[i,])
}
ENA_ts <- ts(Yts, start = (1931+39), frequency = 12)
plot(ENA_ts)
yparp<-pear(ENA_ts,ic="bic")
yparp$model.orders
yparp$phi
yparp$se.phi
 
# ESTUDAR O PERIODO 5 e 10
 
Y_parp <- pear(Yy, ic = "bic")
 
 
Yts <- ts(YY,start=1931,frequency=12)
Xts <- ts(XX,start=1931,frequency=12)
 
if(mes <= 11) {y<-YtM[ ,(mes+1)]; x<-YtM[,mes];nmes<-mes+1;}
 
Y <- as.matrix(dado[5])
typeof(Y)
X <- as.matrix(dado[2:3])
linreg <- lm(Y~X)
summary(linreg)
 
XX <- as.matrix(dado[2])
 
ZZ <- as.matrix(cbind(XX, dado[3]))
 
plot(ENA_ts)
plot(ENA_ts[(40*12+1):length(ENA_ts)], type = "l")
 
ENApad <- scale(ENA, scale = FALSE)
 
 
 
colMeans(ENApad)
apply(ENApad, 2, sd, na.rm = TRUE)
apply(ENApad, 2, mean, na.rm = TRUE)
colMeans(ENApad)
apply(ENA, 2, sd, na.rm = TRUE)
 
 
saida<-"output.txt";
imsaida<-"imoutput.txt";
graph.out<-"./Figures/";
 
# abre arquivo de saida
cabec <- "posto  Nome           t-test var-test Wilcoxon Spearman       MK      ADF       PP    ADF_a     PP_a        mean1        mean2          sd1          sd2";
saida <- paste(arquivo,".out",sep="");
write(cabec,file=saida);
#
ENA<-as.matrix(dado[,2:13])
totAnos <- length(ENA[,1]); # comprimento total da amostra
 
 
 
# ====================================================================
# Xt = ENA como objeto ts a partir da matriz de dados
# Yt = ENA dessazonalizada pela media mensal
# ====================================================================
 
 
((30 - 1) %% 12) + 1
((30 - 1) %/% 12) + 1
 
identical(Yts, tmp)
all.equal(Yts, tmp)
identical(ENA_ts, tmp_ts, num.eq = TRUE)
parp1<-pear(tmp_ts,ic="bic")
parp1$model.orders
parp1$phi
 
 
 
# ====================================================================
# ====================================================================
 
library(partsm)
 
detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0)
ordMax <- 4
aic <- bic <- Fnextp <- Fpval <- rep(NA, ordMax)
lgergnp <- ENA_ts
for(p in 1:ordMax){
  lmpar <- fit.ar.par(wts=lgergnp, detcomp=detcomp, type="PAR", p=p)
  aic[p] <- AIC(lmpar@lm.par, k=2)
  bic[p] <- AIC(lmpar@lm.par, k=log(length(residuals(lmpar@lm.par))))
  Fout <- Fnextp.test(wts=lgergnp, detcomp=detcomp, p=p, type="PAR")
  Fnextp[p] <- Fout@Fstat
  Fpval[p] <- Fout@pval
}
 
aic
bic
Fnextp
Fpval
 
detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0)
ordMax <- 4
aic <- bic <- Fnextp <- Fpval <- rep(NA, ordMax)
#data("gergnp")
#lgergnp <- log(gergnp, base=exp(1))
lgergnp <- ENA_ts
for(p in 1:ordMax){
  lmpar <- fit.ar.par(wts=lgergnp, detcomp=detcomp, type="PAR", p=p)
  aic[p] <- AIC(lmpar@lm.par, k=2)
  bic[p] <- AIC(lmpar@lm.par, k=log(length(residuals(lmpar@lm.par))))
  Fout <- Fnextp.test(wts=lgergnp, detcomp=detcomp, p=p, type="PAR")
  Fnextp[p] <- Fout@Fstat
  Fpval[p] <- Fout@pval
}
 
attr(ENA_ts, "tsp")
out <- pepacf(ENA_ts, 15, plot=TRUE)
 
# =======================================================================================================
# cat(paste("\n"), file=saida, append=TRUE)
# cat(formatC(month.abb,format="s",width=6), file=saida, append=TRUE)
# cat(paste("\n"), file=saida, append=TRUE)
# cat(paste(formatC(parp$model.orders,digits=0,format="f",width=6)), file=saida, append=TRUE)
# cat(paste("\n"))
# 
# cat(parp$model.orders,file = saida, append = TRUE)
 
#cat(title,file=saida,append=T);
#cat(" month",paste(formatC("lag_",format="s",width=11),formatC(1:20,format="d",flag="0",width=2),sep=""),file=saida,fill=13*22+6,append=T);
# cat(formatC(month.abb,format="s",width=6),file=saida,append=TRUE)
# cat("\n",file=saida,append=TRUE)
# cat(formatC(outpepacf$maice,digits=0,format="f",width=6),file=saida,append=TRUE)
# cat("\n",file=saida,append=TRUE)
# cat(formatC(outpepacf$mbice,digits=0,format="f",width=6),file=saida,append=TRUE)
# cat("\n",file=saida,append=TRUE)
# #cat(formatC(month.abb,format="s",width=6));
# #cat(formatC(out$maice,digits=0,format="f",width=6));
# #cat(formatC(out$mbice,digits=0,format="f",width=6))
#
# cat(parp$model.orders, file = saida, append = T)
 
# paste(string_var, collapse = "")
 
# A <- matrix(c(outpepacf$maice,outpepacf$mbice),nrow=2,
#             dimnames = list(c("AIC","BIC"),month.abb))
# write.table(x = A,file = saida,append = T, row.names = T,col.names = T)
# write.table(outpepacf$maice,file = saida,append = T, row.names = T,col.names = T)
 
 
# out <- paste(formatC(month.abb,format="s",width=6),
#              formatC(peacfQt$means,digits=2,format="f",width=13),
#              formatC(peacfQt$standard.deviations,digits=2,format="f",width=13),sep="");
# write(out,file=saida,append=TRUE);