Advertisement
joari

Parp-pear-CCEE

Dec 9th, 2018
166
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 12.05 KB | None | 0 0
  1. # =============================================================================
  2. # r-script para testar ajuste de modelos PAR(p)
  3. # input: Energias medias mensais historicas para REEs (ENA_NE.csv, ...)
  4. # output: relatorio com ordens PAR(p) escolhida via AIC e BIC
  5. # Ref: Estacionaridade das energias afluentes aos REEs
  6. # =============================================================================
  7. # Remove all objects from the current workspace:
  8. rm(list=ls());
  9. #
  10. # =============================================================================
  11. # Load libraries
  12. library(pear)
  13. #
  14. # =============================================================================
  15. # Functions for this study
  16. ajusta_parp <- function(ENA, anoFim, saida){
  17.   tmp <- as.vector(t(ENA))
  18.   ENA_ts <- ts(tmp, end = anoFim, frequency = 12)
  19.   plot(ENA_ts, xlim = c(1931,2018))
  20.   #summary(ENA_ts)
  21.   #hist(ENA_ts,prob=TRUE)
  22.   #boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
  23.  
  24.   #outpepacf <- pepacf(ENA_ts, 15, plot=TRUE)
  25.   parpIC <- list()
  26. #  for (criterion in c("aic", "bic")){
  27.   for (criterion in c("bic")){
  28.     # parp<-pear(ENA_ts,ic=criterion)
  29.     parp<-pear(ENA_ts,ic="bic")
  30.     parpIC <- c(parpIC , parp)
  31.     #parp$model.orders
  32.     IC <- toupper(criterion)
  33.     #title <- paste("\n", IC," model order\n", sep = "" )
  34.     #cat(title, file=saida, append=TRUE)
  35.     #cat(formatC(month.abb,format="s",width=6), file=saida, append=TRUE)
  36.     #cat(paste("\n"), file=saida, append=TRUE)
  37.     cat(paste(formatC(parp$model.orders,digits=NULL,format="d",width=6), sep=""), file=saida, append=TRUE)
  38.     # title <- paste("\n", IC," phi coefficients\n", sep = "" )
  39.     # prt.phi(title, parp$phi, parp$model.orders, saida)  
  40.     # title <- paste(IC," phi standard deviation \n", sep = "" )
  41.     # prt.phi(title, parp$se.phi, parp$model.orders, saida)
  42.    
  43.   }
  44. return(parpIC)
  45. }
  46. # =============================================================================
  47. # print PAR(p) phi or sd.phi coefficients as formatted table
  48. prt.phi <- function(title,phi, model.orders, saida){
  49.   cat(title,file=saida,append=T);
  50.   ordMax <- max(model.orders)
  51.   cat(" month", paste(formatC("lag_",format="s",width=11), formatC(1:ordMax,format="d",flag="0",width=2),sep=""),
  52.       file=saida, fill=13*22+6, append=T);
  53.   for (i in 1:12){
  54.     cat(formatC(phi[i,1:model.orders[[i]]],digits=9,format="f",width=13), labels=paste(formatC(month.abb[i],format="s",width=6)),
  55.         file=saida, fill=13*22+6, append=T);
  56.   }
  57. }
  58. # =============================================================================
  59. #
  60. # ##############################################################################
  61. # Clean workspace & Set working directory
  62. rm(list=ls());                        # Remove all objects from the current workspace:
  63. # workdir <- "D:/ProjetosONS/Cursos/Alex/R-project/ENAs/2018/"; # diretorio de trabalho
  64. workdir <- "~/ProjetosONS/Cursos/Alex/R-project/ENAs/2018/"; # diretorio de trabalho
  65. setwd(workdir);
  66. #
  67. # ==================================================================================
  68. # le dados de ENAs historicas medias mensais
  69. arquivo <- "ENA_NE_2016"
  70. # arqin <- paste(arquivo, ".csv", sep = "")
  71. saida <- paste(arquivo, ".out", sep = "")
  72.  
  73. # dado <- read.csv(arqin, header = T, sep = ";") # tabcsv
  74. #
  75. # ==================================================================================
  76.  
  77. # fit parp (pear package) library("pear")
  78. # ENA <- as.matrix(dado[,2:13])
  79. ENA <- Btt[7,,2:13]
  80. objtmp <- list()
  81. # totAnos <- length(ENA[,1]); # comprimento total da amostra
  82. # anoFim <- 2018
  83. cat("\n BIC model order\n", file=saida, append=TRUE)
  84. title <- paste(" Anoini Anofin   Total", paste(formatC(month.abb,format="s",width=6), collapse = " "), sep = "")
  85. cat(title, file=saida, append=TRUE)
  86. for (inicio in seq(1, 51, by = 5)){
  87. # for (inicio in seq(1, 41, by = 40)){
  88.   count <- count + 1
  89.   anoini <- anoIni + inicio - 1
  90.   out <- paste("\n")
  91.   out <- paste(out,
  92.                formatC(anoini, digits=0,format="f",width=6),
  93.                formatC(anoFim, digits=0,format="f",width=6),
  94.                formatC(anoFim - anoini + 1, digits=0,format="f",width=7), sep = " ")
  95.   cat(out, file=saida, append=TRUE)
  96.   tmp <- ENA[inicio:totYears,]
  97.   objtmp <- c(objtmp , ajusta_parp (tmp, anoFim, saida))
  98. }
  99.  
  100. aa <- array()
  101. for (i in seq(1,11)) {
  102.   nObj <- (i-1)*9  + 1
  103.   print(objtmp[[nObj]][[2]])
  104.   aa[i] <- objtmp[[nObj]][[2]]
  105. }
  106.  
  107. print(aa)
  108. # guradar valores do total de anos para x, ou inicio da serie considerada
  109. plot(x = (anoIni + seq(1, 51, by = 5) - 1), y=aa, main=month.abb[2], ylab = "order")
  110. barplot(aa, names.arg = (anoIni + seq(1, 51, by = 5) - 1),
  111.         xlab = "ano inicial", ylab = "model order", main = paste(arquivo, month.abb[2]))
  112.  
  113.  
  114. seq(1, 41, by = 40)
  115. month.abb[2]
  116. ENA <- as.matrix(dado[,2:13])
  117. totAnos <- length(ENA[,1])
  118.  
  119. tmp <- as.vector(t(ENA))
  120. #tmp <- as.vector(t(ENA[20:85, 1:12]))
  121.  
  122. ENA_ts <- ts(tmp, end = 2018, frequency = 12)
  123.  
  124.  
  125. parp<-pear(ENA_ts,ic="aic")
  126. model <-data.frame()
  127. model <- parp
  128. mydata <- data.frame(parp)
  129. model[1] <- parp
  130. parp$model.orders
  131. parp$phi
  132. parp$se.phi
  133.  
  134. out <- pepacf(ENA_ts, 15, plot=TRUE)
  135. pepacf(ENA_ts, 15, plot=TRUE)
  136. plot(ENA_ts)
  137. summary(ENA_ts)
  138. hist(ENA_ts,prob=TRUE)
  139. boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
  140. acf(as.numeric(ENA_ts))
  141. pacf(as.numeric(ENA_ts))
  142.  
  143.  
  144.  
  145.  
  146. inicio <- 1
  147. anoini <- 1931 + inicio - 1
  148. title <- paste("\n Ano inicial: ", anoini, " Ano final: ", anoFim, " Total : ", anoFim - anoini, sep = "" )
  149. cat(title, file=saida, append=TRUE)
  150. tmp <- ENA[inicio:totAnos,]
  151. ajusta_parp (tmp, anoFim, saida)
  152.  
  153. inicio <- 20
  154. anoini <- 1931 + inicio - 1
  155. title <- paste("\n Ano inicial: ", anoini, " Ano final: ", anoFim, " Total : ", anoFim - anoini, sep = "" )
  156. cat(title, file=saida, append=TRUE)
  157. ENA <- as.matrix(dado[20:totAnos,2:13])
  158. ajusta_parp (ENA, 2018, saida)
  159.  
  160. totAnos <- length(ENA[,1]); # comprimento total da amostra
  161.  
  162. anoFim <- 2016
  163.  
  164. # ==================================================================================
  165. #
  166. # fit parp (pear package) library("pear")
  167. # library(pear)
  168. # ENA_ts = ENA como objeto ts a partir da matriz de dados
  169. # Yts <- ENA[1,]
  170. # for (i in 2:totAnos){
  171. #   Yts <- c(Yts, ENA[i,])
  172. # }
  173.  
  174. tmp <- as.vector(t(ENA))
  175.  
  176. #ENA_ts <- ts(Yts, start = 1931, frequency = 12)
  177. ENA_ts <- ts(tmp, start = 1931, frequency = 12)
  178.  
  179. plot(ENA_ts)
  180. summary(ENA_ts)
  181. hist(ENA_ts,prob=TRUE)
  182. boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
  183. acf(as.numeric(ENA_ts))
  184. pacf(as.numeric(ENA_ts))
  185.  
  186. parp<-pear(ENA_ts,ic="aic")
  187. parp$model.orders
  188. parp$phi
  189. parp$se.phi
  190.  
  191. parp<-pear(ENA_ts, ic="bic")
  192. parp$model.orders
  193.  
  194. # ==================================================================================
  195.  
  196. arquivo <- "ENA_NE_PMO2017"
  197. arqin <- paste(arquivo, ".csv", sep = "")
  198. dado <- read.csv(arqin, header = T, sep = ";") # tabcsv
  199. #
  200. # ==================================================================================
  201.  
  202. ENA <- as.matrix(dado[,2:13])
  203. totAnos <- length(ENA[,1])
  204.  
  205. tmp <- as.vector(t(ENA))
  206. tmp <- as.vector(t(ENA[20:85, 1:12]))
  207.  
  208. ENA_ts <- ts(tmp, end = 2015, frequency = 12)
  209.  
  210. out <- pepacf(ENA_ts, 15, plot=TRUE)
  211. pepacf(ENA_ts, 15, plot=TRUE)
  212. plot(ENA_ts)
  213. summary(ENA_ts)
  214. hist(ENA_ts,prob=TRUE)
  215. boxplot(ENA_ts ~ cycle(ENA_ts), names=month.abb)
  216. acf(as.numeric(ENA_ts))
  217. pacf(as.numeric(ENA_ts))
  218.  
  219. parp<-pear(ENA_ts,ic="aic")
  220. parp$model.orders
  221. parp$phi
  222. parp$se.phi
  223.  
  224. parp<-pear(ENA_ts, ic="bic")
  225. parp$model.orders
  226.  
  227.  
  228.  
  229.  
  230.  
  231. all.equal(parp$model.phi, parp1$model.phi)
  232.  
  233. all(Yts==tmp)
  234. all.equal(Yts,tmp, check.attributes=FALSE)
  235.  
  236. Yts <- ENA[40,]
  237. for (i in 41:totAnos){
  238.   Yts <- c(Yts, ENA[i,])
  239. }
  240. ENA_ts <- ts(Yts, start = (1931+39), frequency = 12)
  241. plot(ENA_ts)
  242. yparp<-pear(ENA_ts,ic="bic")
  243. yparp$model.orders
  244. yparp$phi
  245. yparp$se.phi
  246.  
  247. # ESTUDAR O PERIODO 5 e 10
  248.  
  249. Y_parp <- pear(Yy, ic = "bic")
  250.  
  251.  
  252. Yts <- ts(YY,start=1931,frequency=12)
  253. Xts <- ts(XX,start=1931,frequency=12)
  254.  
  255. if(mes <= 11) {y<-YtM[ ,(mes+1)]; x<-YtM[,mes];nmes<-mes+1;}
  256.  
  257. Y <- as.matrix(dado[5])
  258. typeof(Y)
  259. X <- as.matrix(dado[2:3])
  260. linreg <- lm(Y~X)
  261. summary(linreg)
  262.  
  263. XX <- as.matrix(dado[2])
  264.  
  265. ZZ <- as.matrix(cbind(XX, dado[3]))
  266.  
  267. plot(ENA_ts)
  268. plot(ENA_ts[(40*12+1):length(ENA_ts)], type = "l")
  269.  
  270. ENApad <- scale(ENA, scale = FALSE)
  271.  
  272.  
  273.  
  274. colMeans(ENApad)
  275. apply(ENApad, 2, sd, na.rm = TRUE)
  276. apply(ENApad, 2, mean, na.rm = TRUE)
  277. colMeans(ENApad)
  278. apply(ENA, 2, sd, na.rm = TRUE)
  279.  
  280.  
  281. saida<-"output.txt";
  282. imsaida<-"imoutput.txt";
  283. graph.out<-"./Figures/";
  284.  
  285. # abre arquivo de saida
  286. cabec <- "posto  Nome           t-test var-test Wilcoxon Spearman       MK      ADF       PP    ADF_a     PP_a        mean1        mean2          sd1          sd2";
  287. saida <- paste(arquivo,".out",sep="");
  288. write(cabec,file=saida);
  289. #
  290. ENA<-as.matrix(dado[,2:13])
  291. totAnos <- length(ENA[,1]); # comprimento total da amostra
  292.  
  293.  
  294.  
  295. # ====================================================================
  296. # Xt = ENA como objeto ts a partir da matriz de dados
  297. # Yt = ENA dessazonalizada pela media mensal
  298. # ====================================================================
  299.  
  300.  
  301. ((30 - 1) %% 12) + 1
  302. ((30 - 1) %/% 12) + 1
  303.  
  304. identical(Yts, tmp)
  305. all.equal(Yts, tmp)
  306. identical(ENA_ts, tmp_ts, num.eq = TRUE)
  307. parp1<-pear(tmp_ts,ic="bic")
  308. parp1$model.orders
  309. parp1$phi
  310.  
  311.  
  312.  
  313. # ====================================================================
  314. # ====================================================================
  315.  
  316. library(partsm)
  317.  
  318. detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0)
  319. ordMax <- 4
  320. aic <- bic <- Fnextp <- Fpval <- rep(NA, ordMax)
  321. lgergnp <- ENA_ts
  322. for(p in 1:ordMax){
  323.   lmpar <- fit.ar.par(wts=lgergnp, detcomp=detcomp, type="PAR", p=p)
  324.   aic[p] <- AIC(lmpar@lm.par, k=2)
  325.   bic[p] <- AIC(lmpar@lm.par, k=log(length(residuals(lmpar@lm.par))))
  326.   Fout <- Fnextp.test(wts=lgergnp, detcomp=detcomp, p=p, type="PAR")
  327.   Fnextp[p] <- Fout@Fstat
  328.   Fpval[p] <- Fout@pval
  329. }
  330.  
  331. aic
  332. bic
  333. Fnextp
  334. Fpval
  335.  
  336. detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0)
  337. ordMax <- 4
  338. aic <- bic <- Fnextp <- Fpval <- rep(NA, ordMax)
  339. #data("gergnp")
  340. #lgergnp <- log(gergnp, base=exp(1))
  341. lgergnp <- ENA_ts
  342. for(p in 1:ordMax){
  343.   lmpar <- fit.ar.par(wts=lgergnp, detcomp=detcomp, type="PAR", p=p)
  344.   aic[p] <- AIC(lmpar@lm.par, k=2)
  345.   bic[p] <- AIC(lmpar@lm.par, k=log(length(residuals(lmpar@lm.par))))
  346.   Fout <- Fnextp.test(wts=lgergnp, detcomp=detcomp, p=p, type="PAR")
  347.   Fnextp[p] <- Fout@Fstat
  348.   Fpval[p] <- Fout@pval
  349. }
  350.  
  351. attr(ENA_ts, "tsp")
  352. out <- pepacf(ENA_ts, 15, plot=TRUE)
  353.  
  354. # =======================================================================================================
  355. # cat(paste("\n"), file=saida, append=TRUE)
  356. # cat(formatC(month.abb,format="s",width=6), file=saida, append=TRUE)
  357. # cat(paste("\n"), file=saida, append=TRUE)
  358. # cat(paste(formatC(parp$model.orders,digits=0,format="f",width=6)), file=saida, append=TRUE)
  359. # cat(paste("\n"))
  360. #
  361. # cat(parp$model.orders,file = saida, append = TRUE)
  362.  
  363. #cat(title,file=saida,append=T);
  364. #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);
  365. # cat(formatC(month.abb,format="s",width=6),file=saida,append=TRUE)
  366. # cat("\n",file=saida,append=TRUE)
  367. # cat(formatC(outpepacf$maice,digits=0,format="f",width=6),file=saida,append=TRUE)
  368. # cat("\n",file=saida,append=TRUE)
  369. # cat(formatC(outpepacf$mbice,digits=0,format="f",width=6),file=saida,append=TRUE)
  370. # cat("\n",file=saida,append=TRUE)
  371. # #cat(formatC(month.abb,format="s",width=6));
  372. # #cat(formatC(out$maice,digits=0,format="f",width=6));
  373. # #cat(formatC(out$mbice,digits=0,format="f",width=6))
  374. #
  375. # cat(parp$model.orders, file = saida, append = T)
  376.  
  377. # paste(string_var, collapse = "")
  378.  
  379. # A <- matrix(c(outpepacf$maice,outpepacf$mbice),nrow=2,
  380. #             dimnames = list(c("AIC","BIC"),month.abb))
  381. # write.table(x = A,file = saida,append = T, row.names = T,col.names = T)
  382. # write.table(outpepacf$maice,file = saida,append = T, row.names = T,col.names = T)
  383.  
  384.  
  385. # out <- paste(formatC(month.abb,format="s",width=6),
  386. #              formatC(peacfQt$means,digits=2,format="f",width=13),
  387. #              formatC(peacfQt$standard.deviations,digits=2,format="f",width=13),sep="");
  388. # write(out,file=saida,append=TRUE);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement