Clusterización, perfiles de riesgo y desempeño

##Paquetes

  # install.packages(c("readxl","tibble","dplyr","tidyr","cluster","psych","stringr")) ###Instalar librerias

  ##Librerías
  library(readxl)
  library(tibble)
  library(dplyr)
  library(tidyr)
  library(cluster)
  library(psych)
  library(stringr)

  ##Cargar datos
  ruta_archivo = "C:\\Users\\usuario\\Desktop\\TFG\\Datos\\Base de datos-Para Panel.xlsx"
  data_total   = read_excel(path = ruta_archivo, sheet = "Data")

  ##Filtros base
  entidades_filtrar = c("BAC SAN JOSE","CATHAY","GRUPO MUTUAL","MUCAP","SCOTIABANK","POPULAR")
  data_total$Periodo = as.Date(data_total$Periodo)

  # Filtra por entidades y año (ajusta si deseas más)
  data_filtrada <- data_total[
    data_total$Entidad %in% entidades_filtrar &
      format(data_total$Periodo, "%Y") %in% c("2024"),
  ]

  ##Indicador  Código (I1–I15) 
  #Nota: Aca se el asigna a cada indicar de la SUGEF un codigo de indicar, para fines de simplicidad 
  data_filtrada <- data_filtrada %>%
    mutate(
      Indicador = str_squish(Indicador),
      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 %in% c("Compromiso patrimonial","Solvencia y Capitalización","Solvencia y Capitalizacion") ~ "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
      )
    )

  ##Indicadores a usar en la clusterizacion (k-means)
  #Nota: estos se selecionaron dado temas de liquidez, mora, riesgo etc.
  indicadores_seleccionados <- c("I8","I14","I3","I10","I15")

  ## Media geométrica por Entidad-Indicador
  #Nota: se aplica una media geometrica de los datos por indicador, se usa esta media por su solidez matematica 
  media_geom <- function(x) {
    x <- suppressWarnings(as.numeric(x))
    x <- x[is.finite(x) & x > 0]    # evita log de no-positivos
    if (length(x) == 0) return(NA_real_)
    exp(mean(log(x)))}

  media_geometrica <- data_filtrada %>%
    group_by(Entidad, Indicador) %>%
    summarise(Media_Geometrica = media_geom(Valor), .groups = "drop") %>%
    filter(Indicador %in% indicadores_seleccionados)

  ##k-means (una fila por Entidad)
  df_wide <- media_geometrica %>%
    pivot_wider(names_from = Indicador, values_from = Media_Geometrica)

  if (!"Entidad" %in% names(df_wide)) {
    stop("No hay columna 'Entidad' después del pivot_wider. Verifica datos/indicadores.")
  }

  # Matriz numérica con Entidad 
  mat <- df_wide %>%
    select(Entidad, all_of(indicadores_seleccionados)) %>%  # asegura solo indicadores deseados
    column_to_rownames("Entidad") %>%
    mutate(across(everything(), as.numeric)) %>%
    as.data.frame()

  if (nrow(mat) < 2) stop("Se necesitan al menos 2 entidades para ejecutar k-means/silhouette.")

  # Imputar NAs por media de columna
  mat[is.na(mat)] <- apply(mat, 2, function(x) mean(x, na.rm = TRUE))[col(mat)][is.na(mat)]

  # Eliminar columnas constantes (sd = 0) que rompen el escalado
  sds <- apply(mat, 2, sd, na.rm = TRUE)
  const_cols <- names(sds)[which(!is.finite(sds) | sds == 0)]
  if (length(const_cols) > 0) {
    message("Eliminando columnas sin varianza: ", paste(const_cols, collapse = ", "))
    mat <- mat[, setdiff(colnames(mat), const_cols), drop = FALSE]
  }

  if (ncol(mat) == 0) stop("No quedan indicadores con varianza para clusterizar.")

  # Estandarizar
  mat_num    <- as.matrix(mat)
  mat_scaled <- scale(mat_num)

  if (!all(is.finite(mat_scaled))) {
    stop("Hay NA/NaN/Inf después de escalar. Revisa datos o columnas sin varianza.")
  }

  # Filas útiles y distintas (para limitar k)
  n_rows <- nrow(mat_scaled)
  n_dist <- nrow(unique(as.data.frame(mat_scaled)))
  if (is.null(n_rows) || n_rows < 1) stop("mat_scaled no tiene filas válidas.")

  ##MÉTODO DEL CODO
  set.seed(123)
  max_k <- min(10, n_rows, n_dist)   # tope seguro
  wss <- sapply(1:max_k, function(k){
    tryCatch(
      kmeans(mat_scaled, centers = k, nstart = 25)$tot.withinss,
      error = function(e) NA_real_
    )
  })
  if (all(is.na(wss))) stop("No fue posible calcular WSS para ningún k. Verifica los datos.")

  plot(1:max_k, wss, type = "b", pch = 19,
      xlab = "Número de clusters (k)", ylab = "WSS",
      main = "Método del codo")

  ## MÉTODO SILHOUETTE 
  if (n_rows >= 3 && max_k >= 2) {
    ks_sil <- 2:max_k
    sil_width <- sapply(ks_sil, function(k){
      tryCatch({
        km <- kmeans(mat_scaled, centers = k, nstart = 25)
        mean(silhouette(km$cluster, dist(mat_scaled))[, 3])
      }, error = function(e) NA_real_)
    })
    if (all(is.na(sil_width))) {
      message("No se pudo calcular silhouette para los k probados.")
    } else {
      plot(ks_sil, sil_width, type = "b", pch = 19,
          xlab = "k", ylab = "Silhouette promedio",
          main = "Método Silhouette")
    }
  } else {
    message("Silhouette requiere al menos 3 entidades y k >= 2.")
  }

  #Nota: se incorporan los 2 metodos para veridicar la cnatidad de clusteres


  ##Otros paquetes
  suppressPackageStartupMessages({
    library(dplyr)
    library(tibble)
    library(psych)   # para usar describeBy
  })

  ## Selección k óptimo 
  set.seed(123)

  usar_codo = TRUE   # si quiero usar silhouette poner falso


  if (usar_codo) {
    k_opt =  3   # valor elegido según el gráfico del codo
  } else if (exists("sil_width") && any(is.finite(sil_width))) {
    k_opt = ks_sil[ which.max(sil_width) ]   # mejor k por silhouette
  } else {
    k_opt = min(3, if (exists("max_k")) max_k else 3)  # fallback seguro
  }

  ## Ajuste final del modelo K-means 
  # mat_scaled: matriz/df numérica estandarizada (filas = Entidad)
  stopifnot(exists("mat_scaled"))

  km_fit <- kmeans(
    x         = mat_scaled,
    centers   = k_opt,
    nstart    = 50,
    iter.max  = 100)

  ## RESULTADO PRINCIPAL: Entidad - Cluster 
  # Resolver nombres de entidad de forma robusta
  entidades <- rownames(mat_scaled)
  if (is.null(entidades) || any(is.na(entidades)) || any(entidades == "")) {
    stopifnot(exists("df_wide"), "Entidad" %in% names(df_wide))
    entidades <- df_wide$Entidad}

  cluster_assignments <- tibble(
    Entidad = entidades,
    Cluster = factor(km_fit$cluster, levels = sort(unique(km_fit$cluster))))

  ## Tabla con indicadores + cluster 
  stopifnot(exists("df_wide"), "Entidad" %in% names(df_wide))

  df_clusters <- df_wide %>%
    select(Entidad, everything()) %>%
    left_join(cluster_assignments, by = "Entidad") %>%
    arrange(Cluster, Entidad)

  print(df_clusters, n = Inf)

  ##Conteo por cluster 
  conteo_clusters = df_clusters %>% count(Cluster, name = "Num_Entidades")
  print(conteo_clusters)

  ## Perfil promedio por cluster 
  cluster_profile <- df_clusters %>%
    select(-Entidad) %>%
    group_by(Cluster) %>%
    summarise(
      across(where(is.numeric), mean, na.rm = TRUE, .names = "avg_{.col}"),
      .groups = "drop"
    )

  print(cluster_profile, n = Inf)

  ## (Opcional) Descriptivos por cluster en datos escalados 
  mat_scaled_m <- if (is.data.frame(mat_scaled)) as.matrix(mat_scaled) else mat_scaled
  desc_clusters <- describeBy(mat_scaled_m, group = km_fit$cluster, mat = TRUE)
  print(head(desc_clusters, 10))