Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # =============================================================================
- # 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);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement