ThesiR:

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])]
#alfa.est <- c(.5,.1,.9)
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)


#alfa.est <- c(0.55,.55)

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],
         LONGITUDE=database$LONGITUDE[-id_eta_first],
         LATITUDE=database$LONGITUDE[-id_eta_first]) %>%
  arrange(Grupos)
head(est,13)

#alfa.est <- c(.9,.9,.9)
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
B=300
estt <-matrix(0,nrow(est),B)
for(b in 1:B){
j=1
esti<-NULL
while(j <= 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(500,alfa.aux) %>% mean(trim=.30),
           H=zm*peso,
           betaa=rbeta(1,H,1),
           y = (1-betaa^(delta))^(phi_ts),
           yqcond = (1-(1-quantil)^(1/(bq_t*H)))^(1/aq_t),
           yqmarg = (1-exp(-(delta/peso)*(-log(1-quantil))^(1/alfa.aux)))^phi_ts,
           #ysim= rkumar(n,aq_t,bq_t*zm*peso) %>% median()
           ) %>%
    select(y)

  esti=rbind(esti,e)
  j=j+1
}
  estt[,b]<-as.matrix(esti)
}
head(estt,3)
qest <- apply(estt,1,quantile,prob=c(0.025,0.5,0.975))%>%t()
colnames(qest) <- c("q25","q50","q75")

resul.env <- data.frame(Estações=est$Estações,
                        Grupos=est$Grupos,
                        UR=est$UR,
                        t=est$t,
                        qest) %>% mutate(d=ifelse(UR>=q25 & UR<=q75,1,0))
head(resul.env)
resul.env %>% group_by(Estações) %>% summarise(sum(d)/(TT-1))


fig2 = resul.env %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.)+geom_line(aes(x=t,y=UR))+#geom_point(aes(x=t,y=UR,colour=factor(d)))+
      #geom_line(aes(x=t,y=q50),col=2)+
      #geom_line(aes(x=t,y=q50,colour="red"))+
      geom_ribbon(aes(x=t,ymin=q25,ymax=q75),alpha = 0.5)+
      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(legend.position = "none",axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig2

resul.env.g <- resul.env %>%filter(Grupos=="Grupo 1")

f = resul.env.g %>%
ggplot(.)+geom_line(aes(x=t,y=UR))+#geom_point(aes(x=t,y=UR,colour=factor(d)))+
  #geom_line(aes(x=t,y=q50),col=2)+
  #geom_line(aes(x=t,y=q50,colour="red"))+
  geom_ribbon(aes(x=t,ymin=q25,ymax=q75),alpha = 0.5)+
  ylab("Relative humidity")+ #scale_y_continuous(labels=scales::percent)+
  #scale_x_discrete(labels = c(month.abb))+
  facet_wrap(~Estações,scale="free",ncol=2)+
  theme_bw()+
  #geom_hline(data = d,aes(yintercept = kp))+
  #geom_hline(yintercept = q)+
  theme(legend.position = "none",axis.text.x=element_text(angle=90,size=8))
x11()
f


fig2 = resul.env %>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.)+geom_line(aes(x=t,y=UR))+
      geom_line(aes(x=t,y=q50,colour="red"))+#geom_point(aes(x=t,y=UR,colour=factor(d)))+
      #geom_line(aes(x=t,y=q50),col=2)+
      #geom_line(aes(x=t,y=q50,colour="red"))+
      #geom_ribbon(aes(x=t,ymin=q25,ymax=q75),alpha = 0.5)+
      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(legend.position = "none",axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=3)
x11()
fig2