R/boot_kumar.R

In leonardobfn/ThesiR:

boot_kumar <- function(idx){

  path_estimation <- "estimations_predicts/estimations.rds"
  x <- readRDS(path_estimation)
  #cat( "Current index: ", idx, "\n" )
  formula <- x$formula
  data <- x$data
  w <- x$w
  quantil <- x$quantil
  TT <- x$TT

  mf <- model.frame(Formula::Formula(formula),data=data)
  y <- model.response(mf)
  cov_kl <- model.matrix(Formula::Formula(formula), data = data, rhs = 1)
  cov_delta <- model.matrix(Formula::Formula(formula), data = data, rhs = 2)


  ar <-  x$ar$Estimate
  betas <- x$kl$Estimate
  lambdas <- x$delta$Estimate
  alfas <- x$alpha$Estimate

  link_kl <- x$link.kl
  link_delta <- x$link.delta

  id_eta_first <- seq(1,nrow(data),TT)
  id_eta_last <- seq(TT,nrow(data),TT)

  eta_1 <- cov_kl[,-1]%*%betas[-1]
  eta_2 <- ar*(link_kl$linkfun(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)
  delta <- link_delta$linkinv(cov_delta%*%lambdas)

  est <- data.frame(Estações=data$Estações[-id_eta_first],Grupos=data$Grupos[-id_eta_first],
                    klq,
                    delta=delta[-id_eta_first]) %>%
    mutate(bq_t=1/delta,
           phi_ts = log(klq)/log(1-(1-quantil)^(delta)),
           aq_t = 1/phi_ts) %>% arrange(Grupos)



  estt <- NULL

  for(j in 1:length(alfas)){

    G <- paste("Grupo",j)
    w.aux <- w[[j]];L<-length(w.aux)
    alfa.aux <- alfas[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)
    B=100
    estt.b <-matrix(0,nrow(est),B)
    for(b in 1:B){

      esti<-NULL
      for(j in 1:length(alfas)){
        #qest=NULL
        G <- paste("Grupo",j)
        w.aux <- w[[j]];L<-length(w.aux)
        alfa.aux <- alfas[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=.28),
                 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)

      }
      estt.b[,b]<-as.matrix(esti)
    }


qest <- apply(estt.b,1,quantile,prob=c(0.5))
data <- data %>% arrange(Grupos)
#data[-id_eta_first,"UR"] <-estt$y
data[-id_eta_first,"UR"] <-qest
# z1 = x$data %>% filter(Estações=="Manaus") %>% select(UR)
# z2 = data %>%filter(Estações=="Manaus") %>% select(UR)
# plot(1:120,z1$UR[1:120],type="o")
# lines(1:119,z2$UR[1:119],type="o",col=2)
N <- x$N
erro <- x$erro
#tetastart <- c(ar,betas,lambdas,alfas)

mod <- try(gev_kuma(formula=formula,link.kl = link_kl$name,link.delta = link_delta$name,w = w,data = data,N = N,
                grupos = data$Grupos,quantil = quantil,erro = erro),TRUE)

if(class(mod)=="try-error"){
  ar_est <-  rep(NA,length(ar))
  betas_est <- rep(NA,length(betas))
  lambdas_est <- rep(NA,length(lambdas))
  alfas_est <- rep(NA,length(alfas))

}
if(class(mod)!="try-error"){
  ar_est <-  mod$ar
  betas_est <- mod$kl
  lambdas_est <- mod$delta
  alfas_est <- mod$alpha
}

names <- c("ar_1",
          paste0(colnames(cov_kl),"_kl"),paste0(colnames(cov_delta),"_delta"),paste0("alpha_",1:length(alfas)))

estimation <- data.frame(Estimation=unlist(c(ar_est,betas_est,lambdas_est,alfas_est)),Par=names)
path_monte_carlo <- "estimations_predicts/estimation_mc.txt"
write.table(estimation,file = path_monte_carlo,append = T,quote = T,col.names = F)
#return(estimation)

}