# 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);
}
}
# ==================================================================================
#
Comments
0 B
|👍
/👎
0 B
|👍
/👎