Datos de panel de los clientes institucionales del BANHVI

#Paquetes ==================
  install.packages(c("readxl", "dplyr", "plm", "lmtest", "sandwich", "tidyr", "stringr"), dependencies = TRUE)

  library(readxl)
  library(dplyr)
  library(tidyr)
  library(plm)
  library(lmtest)
  library(sandwich)
  library(stringr)

  # ================== Datos ==================
  ruta_archivo <- file.choose()
  df       <- read_excel(path = ruta_archivo, sheet = "Data")
  df_macro <- read_excel(path = ruta_archivo, sheet = "Data2")

  # ================== Limpieza de data ==================
  # 1) Normalizar nombres posibles de columnas macro (TPC, TBP, IPC, IMAE)
  normaliza_nombres <- function(tbl) {
    nm <- names(tbl)
    
    ren <- c(
      TPC   = intersect(nm, c("TPC","TPCoriginal","TPC_original","tpc")),
      TBP   = intersect(nm, c("TBP","TBPoriginal","TBP_original","tbp")),
      IPC   = intersect(nm, c("IPC","IPCoriginal","IPC_original","ipc")),
      IMAE  = intersect(nm, c("IMAE","IMAEoriginal","IMAE_original","imae")))
    
    for (k in names(ren)) {
      if (length(ren[[k]]) >= 1 && !(k %in% nm)) {
        tbl <- dplyr::rename(tbl, !!k := all_of(ren[[k]][1]))}} tbl}

  df_macro <- normaliza_nombres(df_macro)

  # 2) Fecha y llave mensual
  if (!inherits(df$Periodo, "Date"))    df$Periodo    <- as.Date(df$Periodo)
  if (!inherits(df_macro$Periodo, "Date")) df_macro$Periodo <- as.Date(df_macro$Periodo)

  df       <- df %>% mutate(Periodo_Mes = format(Periodo, "%Y-%m"))
  df_macro <- df_macro %>% mutate(Periodo_Mes = format(Periodo, "%Y-%m"))

  # 3) Revision de duplicados por mes, colapsamos por último valor no-NA
  macro_vars <- intersect(names(df_macro), c("TPC","TBP","IPC","IMAE"))
  df_macro_col <- df_macro %>%
    arrange(Periodo) %>%
    group_by(Periodo_Mes) %>%
    summarise(across(all_of(macro_vars), ~ dplyr::last(na.omit(.))), .groups = "drop")

  # 4) Convertir datos en numérico (por si vinieron como texto con símbolos)
  to_num <- function(x) {
    x <- as.character(x)
    x <- gsub("[^0-9,.-]", "", x)   # eliminar símbolos
    x <- gsub(",", ".", x)          # coma -> punto
    suppressWarnings(as.numeric(x))
  }
  df_macro_col <- df_macro_col %>%
    mutate(across(all_of(macro_vars), to_num))

  # Unir bases 
  df <- df %>%
    left_join(df_macro_col, by = "Periodo_Mes")

  # Etiquetar entidades por grupos
  df <- df %>%
    mutate(GRUPO = case_when(
      Entidad %in% c("COOCIQUE","COOPEALIANZA","COOPEANDE","COOPEMEP","COOPENAE","COOPECAJA") ~ "COOPERATIVAS",
      Entidad %in% c("BAC SAN JOSE","CATHAY","GRUPO MUTUAL","MUCAP","SCOTIABANK","POPULAR") ~ "OTROS",
      TRUE ~ NA_character_
    )) %>%
    filter(!is.na(GRUPO))

  # Renombramiento de indicadores financieros
  df_wide <- df %>%
    select(Entidad, Periodo, Indicador, Valor, all_of(macro_vars), GRUPO) %>%
    mutate(Indicador = case_when(
      Indicador == "Activo Productivo / Activo total" ~ "I1",
      Indicador == "Activo Productivo de Intermediación Financiera / Activo Productivo" ~ "I2",
      Indicador == "Activo Productivo de Intermediación Financiera/ Pasivo con Costo  1" ~ "I3",
      Indicador == "Captaciones a plazo con el público / Pasivo con costo" ~ "I4",
      Indicador == "Cartera al día y con atraso de hasta 90 días (excepto cobro judicial) / Activo Productivo" ~ "I5",
      Indicador == "Cartera al día y con atraso hasta 90 días/Cartera total" ~ "I6",
      Indicador == "Compromiso patrimonial" ~ "I7",
      Indicador == "Estimaciones sobre cartera de créditos / Cartera con atraso mayor a 90 días y cobro judicial" ~ "I8",
      Indicador == "Inversiones en títulos valores / Activo Productivo de Intermediación Financiera" ~ "I9",
      Indicador == "Morosidad mayor a 90 días y cobro judicial / Cartera Directa" ~ "I10",
      Indicador == "Obligaciones con entidades financieras del exterior / Pasivo con costo" ~ "I11",
      Indicador == "Obligaciones con entidades financieras del país / Pasivo con costo" ~ "I12",
      Indicador == "Pasivo con costo / Pasivo Total" ~ "I13",
      Indicador == "Rentabilidad nominal sobre Patrimonio Promedio" ~ "I14",
      Indicador == "Utilidad Operacional Bruta / Gastos de Administración  1" ~ "I15",
      TRUE ~ Indicador
    )) %>%
    pivot_wider(names_from = Indicador, values_from = Valor)


  # Paneles por grupo 
  # Asegurar clase Date
  if (!inherits(df_wide$Periodo, "Date")) df_wide$Periodo <- as.Date(df_wide$Periodo)

  # Quitar posibles duplicados, por Entidad y periodo
  df_wide <- df_wide %>% distinct(Entidad, Periodo, .keep_all = TRUE)

  df_coop  <- df_wide %>% filter(GRUPO == "COOPERATIVAS") %>% as.data.frame()
  df_otros <- df_wide %>% filter(GRUPO == "OTROS")         %>% as.data.frame()

  panel_coop  <- pdata.frame(df_coop,  index = c("Entidad","Periodo"))
  panel_otros <- pdata.frame(df_otros, index = c("Entidad","Periodo"))

  # Modelos (RE) con IMAE 
  cat("========= MODELOS PARA COOPERATIVAS (EFECTOS ALEATORIOS) =========\n")
  modelo_coop_I6  <- plm(I6  ~ TPC + TBP + IPC + IMAE, data = panel_coop,  model = "random", na.action = na.omit)
  modelo_coop_I14 <- plm(I14 ~ TPC + TBP + IPC + IMAE, data = panel_coop,  model = "random", na.action = na.omit)

  cat("\n========= MODELOS PARA OTROS (EFECTOS ALEATORIOS) =========\n")
  modelo_otros_I6  <- plm(I6  ~ TPC + TBP + IPC + IMAE, data = panel_otros, model = "random", na.action = na.omit)
  modelo_otros_I14 <- plm(I14 ~ TPC + TBP + IPC + IMAE, data = panel_otros, model = "random", na.action = na.omit)

  # Resúmen por indicadores
  cat("\n\n=========== RESÚMENES COMPLETOS DE LOS MODELOS ==========\n")

  modelos_coop  <- list(I6 = modelo_coop_I6,  I14 = modelo_coop_I14)
  modelos_otros <- list(I6 = modelo_otros_I6, I14 = modelo_otros_I14)

  cat("\n--- MODELOS COOPERATIVAS ---\n")
  for (nombre in names(modelos_coop)) {
    cat(paste0("\n\n>>> Modelo COOP: ", nombre, "\n"))
    print(summary(modelos_coop[[nombre]]))}

  cat("\n--- MODELOS OTROS ---\n")
  for (nombre in names(modelos_otros)) {
    cat(paste0("\n\n>>> Modelo OTROS: ", nombre, "\n"))
    print(summary(modelos_otros[[nombre]]))}



  ##============================================================================================
  # EFECTOS FIJOS - WITHIN
  library(plm)
  library(lmtest)
  library(sandwich)
  library(dplyr)

  # Reconstruir paneles desde df_wide 
  if (!inherits(df_wide$Periodo, "Date")) df_wide$Periodo <- as.Date(df_wide$Periodo)

  df_coop  <- df_wide %>%
    filter(GRUPO == "COOPERATIVAS") %>%
    distinct(Entidad, Periodo, .keep_all = TRUE) %>%
    as.data.frame()

  df_otros <- df_wide %>%
    filter(GRUPO == "OTROS") %>%
    distinct(Entidad, Periodo, .keep_all = TRUE) %>%
    as.data.frame()

  panel_coop  <- pdata.frame(df_coop,  index = c("Entidad", "Periodo"))
  panel_otros <- pdata.frame(df_otros, index = c("Entidad", "Periodo"))

  cat("========= EFECTOS FIJOS (INDIVIDUAL) - COOPERATIVAS =========\n")
  fe_coop_I6  <- plm(I6  ~ TPC + TBP + IPC + IMAE, data = panel_coop,
                    model = "within", effect = "individual", na.action = na.omit)
  fe_coop_I14 <- plm(I14 ~ TPC + TBP + IPC + IMAE, data = panel_coop,
                    model = "within", effect = "individual", na.action = na.omit)

  cat("\n========= EFECTOS FIJOS (INDIVIDUAL) - OTROS =========\n")
  fe_otros_I6  <- plm(I6  ~ TPC + TBP + IPC + IMAE, data = panel_otros,
                      model = "within", effect = "individual", na.action = na.omit)
  fe_otros_I14 <- plm(I14 ~ TPC + TBP + IPC + IMAE, data = panel_otros,
                      model = "within", effect = "individual", na.action = na.omit)

  # Resultados con errores estándar robustos (cluster por entidad)
  cat("\n--- Resultados COOPERATIVAS (robustos) ---\n")
  print(coeftest(fe_coop_I6,  vcov. = vcovHC(fe_coop_I6,  type = "HC1", cluster = "group")))
  print(coeftest(fe_coop_I14, vcov. = vcovHC(fe_coop_I14, type = "HC1", cluster = "group")))

  cat("\n--- Resultados OTROS (robustos) ---\n")
  print(coeftest(fe_otros_I6,  vcov. = vcovHC(fe_otros_I6,  type = "HC1", cluster = "group")))
  print(coeftest(fe_otros_I14, vcov. = vcovHC(fe_otros_I14, type = "HC1", cluster = "group")))






  #######======================================================================================================
  # HAUSMAN: FE vs RE Nota: esto es para evaluar que modelo Efectos fijos o aleatorios representa mejor el modelo
  library(plm)

  dictamen_haus <- function(test) {
    if (is.null(test$p.value) || is.na(test$p.value)) return("p-valor no disponible")
    if (test$p.value < 0.05) "Prefiere FE (RE inconsistente)" else "Prefiere RE (exogeneidad plausible)"}


  #H₀ (nula): RE es consistente (los efectos no observados de cada entidad no están correlacionados con los regresores). Si H₀ es cierta, RE es además más eficiente que FE.
  #H₁ (alternativa): uno de los modelos es inconsistente (típicamente, RE), por correlación entre efectos no observados y regresores → preferir FE.
  cat("\n================ HAUSMAN FE vs RE (phtest) ================\n")

  cat("\n--- COOPERATIVAS: I6 ---\n")
  haus_coop_I6 <- phtest(fe_coop_I6,  modelo_coop_I6)
  print(haus_coop_I6);  cat("Dictamen: ", dictamen_haus(haus_coop_I6),  "\n")

  cat("\n--- COOPERATIVAS: I14 ---\n")
  haus_coop_I14 <- phtest(fe_coop_I14, modelo_coop_I14)
  print(haus_coop_I14); cat("Dictamen: ", dictamen_haus(haus_coop_I14), "\n")

  cat("\n--- OTROS: I6 ---\n")
  haus_otros_I6 <- phtest(fe_otros_I6,  modelo_otros_I6)
  print(haus_otros_I6);  cat("Dictamen: ", dictamen_haus(haus_otros_I6),  "\n")

  cat("\n--- OTROS: I14 ---\n")
  haus_otros_I14 <- phtest(fe_otros_I14, modelo_otros_I14)
  print(haus_otros_I14); cat("Dictamen: ", dictamen_haus(haus_otros_I14), "\n")


  ###############==============================================================================================
  #REVISION ROBUSTA PARAVER SIGNIFICANCIA DE PANEL 
  # Hay o no efecto panel? (RE vs pooled)
  plmtest(modelo_coop_I6,  type="honda") 
  plmtest(modelo_coop_I14, type="honda")
  plmtest(modelo_otros_I6, type="honda")
  plmtest(modelo_otros_I14,type="honda")
  #Nota: explicacion Es un test LM (Honda) para saber si hay efectos de panel (individuales y/o de tiempo) o no.
  #H0: no hay efectos pooled OLS basta.
  #Si p < 0.05: hay efectos usa modelo de panel (RE o FE).

  #REVISION PARA VER SIGNIFICANCIA DE PANEL 

  pool_coop_I6  <- plm(I6  ~ TPC + TBP + IPC + IMAE, data=panel_coop,  model="pooling")
  pFtest(fe_coop_I6, pool_coop_I6)
  # Nota (pFtest FE vs pooled):
  # H0: sin efectos por entidad (pooled suficiente).
  # Resultado: F = 15.07 (df1 = 5, df2 = 572), p = 6.446e-14  -> Rechazo H0.
  # Conclusión: sí hay heterogeneidad entre entidades; pooled OLS no procede.
  #Si el resultado de esta prueba el Hausman es de p = 1, se usa RE.



  #====================================================================================================
  #Errores robustos
  library(lmtest); library(sandwich)
  coeftest(modelo_coop_I6,  vcov=function(x) vcovHC(x, method="arellano", type="HC1", cluster="group"))


  ##============================================================================================
  #Modelos con rezagos 
  library(plm)
  library(lmtest); library(sandwich)

  # Función corta para mostrar R²
  print_r2 <- function(modelo, etiqueta) {
    s <- summary(modelo)
    cat("\nR-cuadrado —", etiqueta, "\n")
    if (is.null(s$r.squared)) {
      cat("No disponible para este modelo.\n")
    } else {
      print(round(s$r.squared, 4))  # within / between / overall
    }
  }

  #=============================================================================================================================
  # COOPERATIVAS (RE) 
  # I6 con solo lag 2 en cada explicativa
  modelo_coop_I6_lag2 <- plm(
    I6 ~ plm::lag(TPC, 4) + plm::lag(TBP, 6) + plm::lag(IPC, 9) + plm::lag(IMAE, 6),
    data = panel_coop, model = "random", na.action = na.omit
  )
  print_r2(modelo_coop_I6_lag2, "COOP | I6 | lag2")
  coeftest(modelo_coop_I6_lag2, vcov = function(x) vcovHC(x, method="arellano", type="HC1", cluster="group"))

  # I14 con solo lag 2 en cada explicativa
  modelo_coop_I14_lag2 <- plm(
    I14 ~ plm::lag(TPC, 4) + plm::lag(TBP, 6) + plm::lag(IPC, 9) + plm::lag(IMAE, 6),
    data = panel_coop, model = "random", na.action = na.omit
  )
  print_r2(modelo_coop_I14_lag2, "COOP | I14 | lag2")
  coeftest(modelo_coop_I14_lag2, vcov = function(x) vcovHC(x, method="arellano", type="HC1", cluster="group"))



  #=============================================================================================================================
  #  OTROS (RE)
  # I6 con solo lag 2 en cada explicativa
  modelo_otros_I6_lag2 <- plm(
    I6 ~ plm::lag(TPC, 2) + plm::lag(TBP, 2) + plm::lag(IPC, 2) + plm::lag(IMAE, 2),
    data = panel_otros, model = "random", na.action = na.omit
  )
  print_r2(modelo_otros_I6_lag2, "OTROS | I6 | lag2")
  coeftest(modelo_otros_I6_lag2, vcov = function(x) vcovHC(x, method="arellano", type="HC1", cluster="group"))

  # I14 con solo lag 2 en cada explicativa
  modelo_otros_I14_lag2 <- plm(
    I14 ~ plm::lag(TPC, 2) + plm::lag(TBP, 2) + plm::lag(IPC, 2) + plm::lag(IMAE, 2),
    data = panel_otros, model = "random", na.action = na.omit
  )
  print_r2(modelo_otros_I14_lag2, "OTROS | I14 | lag2")
  coeftest(modelo_otros_I14_lag2, vcov = function(x) vcovHC(x, method="arellano", type="HC1", cluster="group"))




  ##============================================================================================
  # PRUEBAS DE LOS MODELOS

  #####Heterocedasticidad (Prueba de Breusch-Pagan modificada)
  bptest(modelo_coop_I6)
  bptest(modelo_otros_I6)
  #Autocorrelación serial (Prueba de Wooldridge)
  pwartest(modelo_coop_I6)   # Wooldridge test
  pwartest(modelo_otros_I6)
  # Errores robustos estándar
  coeftest(modelo_otros_I6, vcov = vcovHC(modelo_otros_I6, type = "HC1"))
  # Errores Driscoll-Kraay (si hay autocorrelación en panel)
  coeftest(modelo_otros_I6, vcov = vcovSCC(modelo_otros_I6))
  #Multicolinealidad (VIF
  # Modelo MCO para grupo "cooperativas" y cálculo de VIF
  ols_model_coop <- lm(I6 ~ TPC + TBP + IPC, data = filter(panel_data, GRUPO == "COOPERATIVAS"))
  vif(ols_model_coop)
  # Modelo MCO para grupo "OTROS bancos " y cálculo de VIF
  ols_model_otros <- lm(I6 ~ TPC + TBP + IPC, data = filter(panel_data, GRUPO == "OTROS"))
  vif(ols_model_otros)