ThesiR:

rm(list=ls())

# Packages -------
pack <- c("tidyverse","extraDistr","devtools","Formula","tictoc","betareg","cowplot")

package.check <- lapply(
  pack,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

devtools::load_all() # loading my functions

# Database------

## Loading database------

data("data_1")

## Parse database------
mes = factor(month.abb,levels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))
database = data_1 %>% filter(Estação=="1") %>%
  mutate(Mês = rep(mes,21) %>% rep(14),
         data=paste0(Ano,"/",Mês),
         t_month =seq(1,252) %>% rep(14),
         t_ = seq(1,12) %>% rep(21) %>% rep(14),
         Grupos = rep(paste0("Grupo ",c(1,2,3,3,2,1,1,3,2,3,3,3,1,1)),each=252),
         cost = cos(2*pi*tt/12),sent=sin(2*pi*tt/12))



head(database)


## Calculating the weights ---------

u1 = "Nova Olinda do Norte";u2="Maraã";u3="Itamarati"
u = c(u1,u2,u3)

w = list()
for(j in 1:3){
  g = paste0("Grupo ",j)
  dados_aux = database %>% dplyr::filter(Grupos==g)
  u_lat_log = dplyr::filter(data_1,Estações==u[j])[c("Estações","LATITUDE","LONGITUDE")]
  lat_log = unique(dados_aux[c("Estações","LATITUDE","LONGITUDE")])
  aux = rbind(lat_log,u_lat_log)
  w[j] = weight(mun = aux$Estações,u_m=u[j],lat=as.numeric(aux$LATITUDE),long=as.numeric(aux$LONGITUDE))
}
names(w)=u
w
sapply(w, sum)
head(database)

# Estimation ---------
TT = 120 #
formula <- UR~log(Altitude)+sent+cost|sent+cost
data=database %>% group_by(Estações) %>% slice(1:TT)
nrow(data)
grupos = data$Grupos
w=w
quantil=.5
N=500
erro = 10^(-4)
link.kl="loglog"
link.delta="log"

tictoc::tic()
x <- gev_kuma(formula=formula,
              grupos = data$Grupos,
              w=w,
              N=N,
              quantil=.5,
              data=data,
              link.kl="loglog",link.delta="log",erro=erro
              )
tictoc::toc()
x$u<-u
#1149.24 sec elapsed
path_estimation <- "estimations_predicts/estimations.rds"
saveRDS(x,path_estimation) # saving the results

# Generating the database/Analysis of results with estimates paramenters-------

require(extraDistr)
require(tidyverse)
rm(list=ls())
path_estimation <- "estimations_predicts/estimations.rds"
results <- readRDS(path_estimation)
database <- results$data
cov_kl <- results$cov_kl
cov_delta <- results$cov_delta
y <- results$data$UR
ncx <- ncol(cov_kl)
ncv <- ncol(cov_delta)
grupos<-results$data$Grupos
w <- results$w
quantil<-results$quantil
grupos= as.factor(grupos);gr = levels((grupos))
DADOS = data.frame(database$Estações,y,cov_kl,cov_delta,grupos)
colnames(DADOS) = c("mun","y",
                    colnames(cov_kl),
                    paste0(colnames(cov_delta),"_delta"),
                    "grupos")

EMVteta2 <- c(results$ar$Estimate,
              results$kl$Estimate,
              results$delta$Estimate,
              results$alpha$Estimate)

betas <- EMVteta2[-1][(1:ncx)]
lambda <- EMVteta2[-1][(ncx+1):(ncx+ncv)]
phi_ar <- EMVteta2[1]
alfa.est <- EMVteta2[-1][(ncx+ncv+1):length(EMVteta2[-1])]

link.kl <- results$link.kl
link.delta <- results$link.delta
TT=results$TT
delta <- link.delta$linkinv(as.matrix(DADOS[c(paste0(colnames(cov_delta),"_delta"))]) %*% lambda)
id_eta_first <- seq(1,nrow(DADOS),TT)
id_eta_last <- seq(TT,nrow(DADOS),TT)
eta_1 <- cov_kl[,-1]%*%betas[-1]
eta_2 <- phi_ar*(link.kl$linkfun(DADOS$y) - eta_1)
eta_kl <- betas[1]+eta_1[-c(id_eta_first)] + eta_2[-c(id_eta_last)]
klq <- link.kl$linkinv(eta_kl)




est <- data.frame(Estações=DADOS$mun[-id_eta_first],Grupos=DADOS$grupos[-id_eta_first],
                  klq,
                  delta=delta[-id_eta_first],
                  UR=database$UR[-id_eta_first]) %>%
  mutate(bq_t=1/delta,
         phi_ts = log(klq)/log(1-(1-quantil)^(delta)),
         aq_t = 1/phi_ts,
         t=database$t[-id_eta_first],
         UR=database$UR[-id_eta_first],
         Mês=database$Mês[-id_eta_first],
         Ano=database$Ano[-id_eta_first],
         LONGITUDE=database$LONGITUDE[-id_eta_first],
         LATITUDE=database$LATITUDE[-id_eta_first]) %>%
       arrange(Grupos)
head(est,13)

estt = NULL

rstable_pos <- function(n,alfa){
  U <- runif(n,0,pi)
  v <- rexp(n,1)
  a_U <- ((sin(alfa*U)/sin(U))^(1/(1-alfa)))*((sin((1-alfa)*U))/sin(alfa*U))
  E <- (a_U/v)^((1-alfa)/alfa)
  return(E)

}

n=500
estt <-NULL

for(j in 1:length(alfa.est)){
  #qest=NULL
  G <- paste("Grupo",j)
  w.aux <- w[[j]];L<-length(w.aux)
  alfa.aux <- alfa.est[j]
  #set.seed(10)
  e <- est %>% dplyr::filter(Grupos==G) %>% mutate(peso=w.aux %>% rep(each=(TT-1)),zm = rstable_pos(1000,alfa.aux) %>% mean(trim=0.28),
                                                   H=zm*peso,
                                                   Espy_z= bq_t*H*beta(1+1/aq_t,bq_t*H),
                                                   Vary_z= bq_t*H*beta(1+2/aq_t,bq_t*H)-(bq_t*H*beta(1+1/aq_t,bq_t*H))^2,
                                                   betaa=rbeta(1,H,1),
                                                   y = (1-betaa^(delta))^(phi_ts),
                                                   ycond = (1-(1-quantil)^(1/(bq_t*H)))^(1/aq_t),
                                                   ymarg = (1-exp(-(delta/peso)*(-log(1-quantil))^(1/alfa.aux)))^phi_ts,
                                                   Fy_z = pkumar(.7,aq_t,bq_t*zm*peso,lower.tail = T),
                                                   cofvar=sqrt(Vary_z)/Espy_z,
                                                  )


  estt=rbind(estt,e)
}
estt <- estt %>% arrange(Grupos)
q_est = NULL
for(i in 1:nrow(estt)){
  #zm <- stabledist::rstable(n=1,alpha=alfa.est[1],beta=1,gamma=1,delta = 0,pm=1)
  #betas <- rbeta(600,e$peso[i]*zm,1)
  # q_est[i] <- (1-betas^(delta[i]))^(e$phi[i]) %>% mean(trim=.10)
  #q_est[i] <- (1-betas^(delta[i]))^(e$phi[i]) %>% mean(trim=.10)
  q_est[i] <- rkumar(n,estt$aq_t[i],estt$bq_t[i]*estt$zm*estt$peso[i]) %>% median()


}

 estt$q_est <- q_est
 head(estt,24)
 #estt %>% group_by(Estações,t) %>% summarise(m=median(q_est))
 est = estt %>% mutate(q_est=y,erro1=round(((UR-q_est)/UR),2),erro2=scale(UR-q_est))
 est

# Figuras para an?lise

fig1 = estt %>% select(t,q_est,ymarg,UR,y,Grupos,Estações)%>%pivot_longer(cols = c(y,UR)) %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.,aes(colour=name  )) + geom_line(aes(t,value))+
      #geom_point(aes(t,value))+
      ylab("Relative humidity")+ #scale_y_continuous(labels=scales::percent)+
      #scale_x_discrete(labels = c(month.abb))+
      facet_grid(~Estações,scale="free")+
      theme_bw()+#geom_hline(data = d,aes(yintercept = kp))+
      #geom_hline(yintercept = q)+
      theme(axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig1

saveRDS(estt,"valores_preditos.rds")

fig2 = est %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.) + geom_boxplot(aes(Mês,erro2))+
      ylab("Relative humidity")+ scale_y_continuous(labels=scales::percent)+
      scale_x_discrete(labels = c(month.abb))+
      facet_grid(~Estações)+geom_hline(yintercept = 0)+
      theme_bw()+
      theme(axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig2



fig3 = est %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.) + geom_point(aes(t,erro2))+
      ylab("Relative humidity")+ scale_y_continuous()+
      # scale_x_discrete(labels = c(month.abb))+
      facet_wrap(~Estações,scale="free_y")+
      theme_bw()+
      theme(axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig3


ic_alpha= function(alpha, acf_res){
  return(qnorm((1 + (1 - alpha))/2)/sqrt(acf_res$n.used))
}

acff <- est %>% group_by(Estações,Grupos) %>% summarise(acff=acf(erro2,plot=F)$acf[,,1],
                                                        lag = acf(erro2,plot=F)$lag,
                                                        lim1=ic_alpha(0.05,acf(erro2,plot=F)))



fig3 = acff %>% data.frame() %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.,aes(x=lag,y=acff)) + geom_hline(aes(yintercept = 0)) +
      geom_segment(mapping = aes(xend = lag, yend = 0))+
      ylab("ACF")+ scale_y_continuous()+
      geom_hline(aes(yintercept = lim1), linetype = 2, color = 'blue') +
      geom_hline(aes(yintercept = -lim1), linetype = 2, color = 'blue')+
      # scale_x_discrete(labels = c(month.abb))+
      facet_grid(~Estações,scale="free_x")+
      theme_bw()+
      theme(axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig3