################################################################################
# PRONOSTICO DE SERIES - MERCADO DE SERVICIOS MOVILES COLOMBIA
# Pradatos - Analitica + IA - Marzo 2026
# Base: ServiciosMoviles.xlsx 
################################################################################

rm(list = ls(all.names = TRUE))

if (interactive()) {
  cat("\n== PRADATOS - Pronostico de Series Moviles ==\n\n")
  ruta_excel <- file.choose()
  ruta_salida_excel <- choose.dir(caption = "Carpeta Excel salida")
  ruta_salida_graficas <- choose.dir(caption = "Carpeta graficas")
} else {
  ruta_excel <- "ServiciosMoviles.xlsx"
  ruta_salida_excel <- "salida"
  ruta_salida_graficas <- "salida/graficas"
}
dir.create(ruta_salida_excel, showWarnings=FALSE, recursive=TRUE)
dir.create(ruta_salida_graficas, showWarnings=FALSE, recursive=TRUE)

paquetes <- c("readxl","compositions","vars","forecast","tseries",
              "ggplot2","dplyr","tidyr","writexl","scales")
for (p in paquetes) { if (!requireNamespace(p,quietly=TRUE)) install.packages(p); library(p,character.only=TRUE) }

options(scipen = 999)

# ═══════════════════════════════════════════════════════════════
# 1. LECTURA Y CONSOLIDACION (identico al codigo original)
# ═══════════════════════════════════════════════════════════════
cat("[1/7] Leyendo y consolidando datos...\n")
Serviciosmoviles <- read_excel(ruta_excel, sheet = 1)

# Consolidar: Movistar+Tigo = "Tigo recargado", Claro, Partners, resto = "Otros"
Serviciosmoviles2 <- Serviciosmoviles %>%
  mutate(EMPRESA = case_when(
    EMPRESA %in% c("Movistar", "Tigo") ~ "Tigo recargado",
    EMPRESA %in% c("Claro", "Partners") ~ EMPRESA,
    TRUE ~ "Otros"))

Serviciosmoviles2 <- Serviciosmoviles2 %>%
  group_by(ANNO, TRIMESTRE, EMPRESA) %>%
  summarise(across(where(is.numeric), sum, na.rm = TRUE), .groups = "drop")

# CRITICO: Partners no existe en periodos tempranos. Completar con 0
all_periodos <- Serviciosmoviles2 %>%
  dplyr::select(ANNO, TRIMESTRE) %>%
  distinct()
all_empresas <- data.frame(EMPRESA = c("Claro","Tigo recargado","Partners","Otros"))
grilla <- tidyr::crossing(all_periodos, all_empresas)
Serviciosmoviles2 <- grilla %>%
  left_join(Serviciosmoviles2, by = c("ANNO","TRIMESTRE","EMPRESA")) %>%
  mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

# 1a. Cuotas de mercado (identico: cada variable / sum por trimestre)
Serviciosmoviles2 <- Serviciosmoviles2 %>%
  group_by(ANNO, TRIMESTRE) %>%
  mutate(across(-any_of(c("ANNO","TRIMESTRE","EMPRESA")),
                ~ .x / sum(.x, na.rm = TRUE), .names = "CUOTA_{.col}")) %>%
  ungroup()

# 1b. Tasas de crecimiento
Serviciosmoviles2 <- Serviciosmoviles2 %>%
  arrange(EMPRESA, ANNO, TRIMESTRE) %>%
  group_by(EMPRESA) %>%
  mutate(across(matches("^(ABONADOS|SUSCRIPTORES|USUARIOS|TRAFICO)"),
                ~ { g <- (.x - lag(.x)) / lag(.x); g[is.infinite(g) | is.nan(g)] <- NA; g },
                .names = "TASAC_{.col}")) %>%
  ungroup()

# 1c. IPROM, ARPU, CONSUMO (identico al codigo original)
Serviciosmoviles2 <- Serviciosmoviles2 %>%
  arrange(EMPRESA, ANNO, TRIMESTRE) %>%
  group_by(EMPRESA) %>%
  mutate(
    IPROM_VM     = INGRESOS_VM / TRAFICO_VM,
    IPROM_VM_PRE = INGRESOS_VM_PRE / TRAFICO_VM_PRE,
    IPROM_VM_POS = INGRESOS_VM_POS / TRAFICO_VM_POS,
    IPROM_IM     = INGRESOS_IM / TRAFICO_IM,
    IPROM_IM_DEM = INGRESOS_IMDEM / TRAFICO_IMDEM,
    IPROM_IM_CF  = INGRESOS_IM_CF / TRAFICO_IM_CF,
    ARPU_VM      = INGRESOS_VM / USUARIOS_VM,
    ARPU_VM_PRE  = INGRESOS_VM_PRE / ABONADOS_VM,
    ARPU_VM_POS  = INGRESOS_VM_POS / SUSCRIPTORES_VM,
    ARPU_IM      = INGRESOS_IM / USUARIOS_IM,
    ARPU_IM_DEM  = INGRESOS_IMDEM / ABONADOS_IM,
    ARPU_IM_CF   = INGRESOS_IM_CF / SUSCRIPTORES_IM,
    ARPU_SM      = INGRESOS_SM / USUARIOS_SM,
    CONS_IM      = TRAFICO_IM / USUARIOS_IM,
    CONS_VM      = TRAFICO_VM / USUARIOS_VM
  ) %>%
  ungroup()

# 1d. Stenbacka (identico al codigo original)
calcular_stenbacka <- function(base) {
  base %>%
    group_by(ANNO, TRIMESTRE) %>%
    mutate(
      S1_USUARIOS_VM = max(CUOTA_USUARIOS_VM, na.rm=TRUE),
      S2_USUARIOS_VM = sort(CUOTA_USUARIOS_VM, decreasing=TRUE)[2],
      STEN_USUARIOS_VM_G1  = 0.5 * (1 - (S1_USUARIOS_VM^2 - S2_USUARIOS_VM^2)),
      STEN_USUARIOS_VM_G05 = 0.5 * (1 - 0.5 * (S1_USUARIOS_VM^2 - S2_USUARIOS_VM^2)),

      S1_TRAFICO_VM = max(CUOTA_TRAFICO_VM, na.rm=TRUE),
      S2_TRAFICO_VM = sort(CUOTA_TRAFICO_VM, decreasing=TRUE)[2],
      STEN_TRAFICO_VM_G1  = 0.5 * (1 - (S1_TRAFICO_VM^2 - S2_TRAFICO_VM^2)),
      STEN_TRAFICO_VM_G05 = 0.5 * (1 - 0.5 * (S1_TRAFICO_VM^2 - S2_TRAFICO_VM^2)),

      S1_INGRESOS_VM = max(CUOTA_INGRESOS_VM, na.rm=TRUE),
      S2_INGRESOS_VM = sort(CUOTA_INGRESOS_VM, decreasing=TRUE)[2],
      STEN_INGRESOS_VM_G1  = 0.5 * (1 - (S1_INGRESOS_VM^2 - S2_INGRESOS_VM^2)),
      STEN_INGRESOS_VM_G05 = 0.5 * (1 - 0.5 * (S1_INGRESOS_VM^2 - S2_INGRESOS_VM^2)),

      S1_USUARIOS_IM = max(CUOTA_USUARIOS_IM, na.rm=TRUE),
      S2_USUARIOS_IM = sort(CUOTA_USUARIOS_IM, decreasing=TRUE)[2],
      STEN_USUARIOS_IM_G1  = 0.5 * (1 - (S1_USUARIOS_IM^2 - S2_USUARIOS_IM^2)),
      STEN_USUARIOS_IM_G05 = 0.5 * (1 - 0.5 * (S1_USUARIOS_IM^2 - S2_USUARIOS_IM^2)),

      S1_TRAFICO_IM = max(CUOTA_TRAFICO_IM, na.rm=TRUE),
      S2_TRAFICO_IM = sort(CUOTA_TRAFICO_IM, decreasing=TRUE)[2],
      STEN_TRAFICO_IM_G1  = 0.5 * (1 - (S1_TRAFICO_IM^2 - S2_TRAFICO_IM^2)),
      STEN_TRAFICO_IM_G05 = 0.5 * (1 - 0.5 * (S1_TRAFICO_IM^2 - S2_TRAFICO_IM^2)),

      S1_INGRESOS_IM = max(CUOTA_INGRESOS_IM, na.rm=TRUE),
      S2_INGRESOS_IM = sort(CUOTA_INGRESOS_IM, decreasing=TRUE)[2],
      STEN_INGRESOS_IM_G1  = 0.5 * (1 - (S1_INGRESOS_IM^2 - S2_INGRESOS_IM^2)),
      STEN_INGRESOS_IM_G05 = 0.5 * (1 - 0.5 * (S1_INGRESOS_IM^2 - S2_INGRESOS_IM^2)),

      S1_USUARIOS_SM = max(CUOTA_USUARIOS_SM, na.rm=TRUE),
      S2_USUARIOS_SM = sort(CUOTA_USUARIOS_SM, decreasing=TRUE)[2],
      STEN_USUARIOS_SM_G1  = 0.5 * (1 - (S1_USUARIOS_SM^2 - S2_USUARIOS_SM^2)),
      STEN_USUARIOS_SM_G05 = 0.5 * (1 - 0.5 * (S1_USUARIOS_SM^2 - S2_USUARIOS_SM^2)),

      S1_INGRESOS_SM = max(CUOTA_INGRESOS_SM, na.rm=TRUE),
      S2_INGRESOS_SM = sort(CUOTA_INGRESOS_SM, decreasing=TRUE)[2],
      STEN_INGRESOS_SM_G1  = 0.5 * (1 - (S1_INGRESOS_SM^2 - S2_INGRESOS_SM^2)),
      STEN_INGRESOS_SM_G05 = 0.5 * (1 - 0.5 * (S1_INGRESOS_SM^2 - S2_INGRESOS_SM^2))
    ) %>% ungroup()
}

Serviciosmoviles2 <- calcular_stenbacka(Serviciosmoviles2)

# 1e. HHI (identico: 10000 * sum(cuota^2))
HHI2 <- Serviciosmoviles2 %>%
  dplyr::group_by(ANNO, TRIMESTRE) %>%
  dplyr::summarise(
    dplyr::across(
      starts_with("CUOTA_"),
      ~ 10000 * sum((.x)^2, na.rm = TRUE),
      .names = "HHI_{.col}"
    ),
    .groups = "drop"
  )

# Agregar periodo
Serviciosmoviles2 <- Serviciosmoviles2 %>%
  mutate(TRIMESTRE_ROM = as.character(as.roman(TRIMESTRE)),
         FECHA = paste0(ANNO, "-", TRIMESTRE_ROM),
         FECHA_NUM = ANNO + (TRIMESTRE - 1) / 4)

HHI2 <- HHI2 %>%
  mutate(TRIMESTRE_ROM = as.character(as.roman(TRIMESTRE)),
         FECHA = paste0(ANNO, "-", TRIMESTRE_ROM),
         FECHA_NUM = ANNO + (TRIMESTRE - 1) / 4)

# Usar esta base como df principal
df <- Serviciosmoviles2
periodos <- df %>% filter(EMPRESA=="Claro") %>% arrange(ANNO,TRIMESTRE) %>% pull(FECHA)
fechas <- df %>% filter(EMPRESA=="Claro") %>% arrange(ANNO,TRIMESTRE) %>% pull(FECHA_NUM)
n_obs <- length(periodos)
empresas <- c("Claro","Tigo recargado","Partners","Otros")

h <- 3
ua <- max(df$ANNO); ut <- max(df$TRIMESTRE[df$ANNO==ua])
ff <- numeric(h); pf <- character(h)
for (i in 1:h) { t<-ut+i; a<-ua; while(t>4){t<-t-4;a<-a+1}; ff[i]<-a+(t-1)/4; pf[i]<-paste0(a,"-",as.character(as.roman(t))) }

cat("  n=",n_obs," h=",h," empresas:", paste(empresas,collapse=", "),"\n")
cat("  Ultimo:",periodos[n_obs]," Hasta:",pf[h],"\n")

# Paleta
PAL <- c("Claro"="#E30613","Tigo recargado"="#0033A0","Partners"="#8B008B","Otros"="#C8B89A","Total"="#1E1E1E")

# Series
cuotas <- c("CUOTA_USUARIOS_VM","CUOTA_TRAFICO_VM","CUOTA_INGRESOS_VM",
  "CUOTA_TRAFICO_VM_PRE","CUOTA_TRAFICO_VM_POS","CUOTA_INGRESOS_VM_PRE","CUOTA_INGRESOS_VM_POS",
  "CUOTA_USUARIOS_IM","CUOTA_TRAFICO_IM","CUOTA_INGRESOS_IM",
  "CUOTA_TRAFICO_IMDEM","CUOTA_TRAFICO_IM_CF","CUOTA_INGRESOS_IMDEM","CUOTA_INGRESOS_IM_CF",
  "CUOTA_ABONADOS_VM","CUOTA_SUSCRIPTORES_VM","CUOTA_USUARIOS_SM","CUOTA_INGRESOS_SM")
volumenes <- c("USUARIOS_VM","TRAFICO_VM","INGRESOS_VM","USUARIOS_IM","TRAFICO_IM","INGRESOS_IM","USUARIOS_SM","INGRESOS_SM")
precios <- c("IPROM_VM","IPROM_IM","ARPU_VM","ARPU_IM","ARPU_SM","CONS_IM","CONS_VM")

ETQ <- c(CUOTA_USUARIOS_VM="Cuota Usuarios Voz",CUOTA_TRAFICO_VM="Cuota Trafico Voz",
  CUOTA_INGRESOS_VM="Cuota Ingresos Voz",CUOTA_TRAFICO_VM_PRE="Cuota Trafico Voz Pre",
  CUOTA_TRAFICO_VM_POS="Cuota Trafico Voz Pos",CUOTA_INGRESOS_VM_PRE="Cuota Ingresos Voz Pre",
  CUOTA_INGRESOS_VM_POS="Cuota Ingresos Voz Pos",CUOTA_USUARIOS_IM="Cuota Usuarios Internet",
  CUOTA_TRAFICO_IM="Cuota Trafico Internet",CUOTA_INGRESOS_IM="Cuota Ingresos Internet",
  CUOTA_TRAFICO_IMDEM="Cuota Trafico Internet Dem",CUOTA_TRAFICO_IM_CF="Cuota Trafico Internet CF",
  CUOTA_INGRESOS_IMDEM="Cuota Ingresos Internet Dem",CUOTA_INGRESOS_IM_CF="Cuota Ingresos Internet CF",
  CUOTA_ABONADOS_VM="Cuota Abonados Voz",CUOTA_SUSCRIPTORES_VM="Cuota Suscriptores Voz",
  CUOTA_USUARIOS_SM="Cuota Usuarios SM",CUOTA_INGRESOS_SM="Cuota Ingresos SM",
  USUARIOS_VM="Usuarios Voz",TRAFICO_VM="Trafico Voz (min)",INGRESOS_VM="Ingresos Voz (COP)",
  USUARIOS_IM="Usuarios Internet",TRAFICO_IM="Trafico Internet (MB)",INGRESOS_IM="Ingresos Internet (COP)",
  USUARIOS_SM="Usuarios Serv Moviles",INGRESOS_SM="Ingresos Serv Moviles",
  IPROM_VM="IPROM Voz ($/min)",IPROM_IM="IPROM Internet ($/MB)",
  ARPU_VM="ARPU Voz",ARPU_IM="ARPU Internet",ARPU_SM="ARPU Servicios Moviles",
  CONS_IM="Consumo Internet (MB/acceso)",CONS_VM="Consumo Voz (min/acceso)")

# ═══════════════════════════════════════════════════════════════
# 2. GRUPO A: ILR + VAR (cuotas)
# ═══════════════════════════════════════════════════════════════
cat("[2/7] Grupo A: Cuotas (ILR+VAR)...\n")
res_cuotas <- list()
for (vc in cuotas) {
  cat("  ->",ETQ[vc],"\n")
  mat <- matrix(NA,n_obs,4); colnames(mat)<-empresas
  for (j in 1:4) {
    v <- df %>% filter(EMPRESA==empresas[j]) %>% arrange(ANNO,TRIMESTRE) %>% pull(!!sym(vc))
    if(length(v)>=n_obs) mat[,j] <- as.numeric(v[1:n_obs])
  }
  mat[is.na(mat)]<-1e-6; mat[mat<=0]<-1e-6; mat<-mat/rowSums(mat)
  cm<-acomp(mat); im<-ilr(cm)
  ml<-max(1,min(4,floor((n_obs-1)/(ncol(im)+1))-1))
  tryCatch({
    sel<-VARselect(im,lag.max=ml,type="const"); po<-max(1,min(sel$selection["SC(n)"],ml))
    mod<-VAR(im,p=po,type="const"); pv<-predict(mod,n.ahead=h,ci=0.95)
    pi<-plo<-phi<-matrix(NA,h,ncol(im))
    for(k in 1:ncol(im)){nm<-colnames(im)[k];if(is.null(nm))nm<-paste0("y",k);pk<-pv$fcst[[nm]];pi[,k]<-pk[,"fcst"];plo[,k]<-pk[,"lower"];phi[,k]<-pk[,"upper"]}
    ps<-as.matrix(ilrInv(pi,orig=cm));sl<-as.matrix(ilrInv(plo,orig=cm));su<-as.matrix(ilrInv(phi,orig=cm))
    colnames(ps)<-colnames(sl)<-colnames(su)<-empresas
    res_cuotas[[vc]]<-list(hist=mat,pred=ps,lo=sl,hi=su,mod=paste0("VAR(",po,")+ILR"))
  },error=function(e){
    cat("    !",e$message,"\n"); pn<-matrix(rep(mat[n_obs,],each=h),nrow=h);colnames(pn)<-empresas
    res_cuotas[[vc]]<<-list(hist=mat,pred=pn,lo=pn*0.9,hi=pmin(pn*1.1,1),mod="Naive")
  })
}

# ═══════════════════════════════════════════════════════════════
# 3. GRUPO B: VOLUMENES (competencia de modelos + CV temporal)
# ═══════════════════════════════════════════════════════════════
cat("[3/7] Grupo B: Volumenes...\n")

test_estacionariedad <- function(y) {
  adf_p <- tryCatch(adf.test(y)$p.value, error=function(e) NA)
  kpss_p <- tryCatch(kpss.test(y, null="Level")$p.value, error=function(e) NA)
  d_sug <- ndiffs(y, test="kpss")
  list(adf_p=adf_p, kpss_p=kpss_p, d=d_sug)
}

cv_temporal <- function(y_ts, h_fc, k_test=6) {
  n <- length(y_ts)
  if (n < 20+h_fc) k_test <- max(2, n-15-h_fc)
  origenes <- seq(n-k_test-h_fc+1, n-h_fc)
  cands <- list(
    ARIMA_auto=function(tr) auto.arima(tr,stepwise=FALSE,approximation=FALSE),
    ETS_auto=function(tr) ets(tr), ETS_damped=function(tr) ets(tr,damped=TRUE),
    ARIMA_drift=function(tr) Arima(tr,order=c(0,1,0),include.drift=TRUE))
  rmse_mat <- matrix(NA,length(origenes),length(cands)); colnames(rmse_mat)<-names(cands)
  for (i in seq_along(origenes)) {
    o<-origenes[i]; train<-window(y_ts,end=c(floor(time(y_ts)[o]),cycle(y_ts)[o]))
    test_r<-y_ts[(o+1):min(o+h_fc,n)]; hr<-length(test_r); if(hr<1)next
    for (m in seq_along(cands)) {
      tryCatch({ fit<-cands[[m]](train); fc<-forecast(fit,h=hr,level=95)
        err<-as.numeric(test_r)-as.numeric(fc$mean[1:hr]); rmse_mat[i,m]<-sqrt(mean(err^2,na.rm=TRUE))
      },error=function(e){})
    }
  }
  rp<-colMeans(rmse_mat,na.rm=TRUE); rp[is.nan(rp)]<-Inf; list(rmse=rp,mejor=names(which.min(rp)))
}

ajustar_ganador <- function(y_ts,nm,hf) {
  if(nm=="ARIMA_auto"){fit<-auto.arima(y_ts,stepwise=FALSE,approximation=FALSE);fc<-forecast(fit,h=hf,level=95);mt<-paste0("ARIMA(",paste(fit$arma[c(1,6,2)],collapse=","),")")
  } else if(nm=="ETS_auto"){fit<-ets(y_ts);fc<-forecast(fit,h=hf,level=95);mt<-fit$method
  } else if(nm=="ETS_damped"){fit<-ets(y_ts,damped=TRUE);fc<-forecast(fit,h=hf,level=95);mt<-fit$method
  } else {fit<-Arima(y_ts,order=c(0,1,0),include.drift=TRUE);fc<-forecast(fit,h=hf,level=95);mt<-"ARIMA(0,1,0)+drift"}
  lb_p<-tryCatch(Box.test(residuals(fc),lag=min(10,length(y_ts)/5),type="Ljung-Box")$p.value,error=function(e)NA)
  list(pred=as.numeric(fc$mean),lo=as.numeric(fc$lower[,1]),hi=as.numeric(fc$upper[,1]),mod=mt,lb_p=lb_p)
}

res_vol <- list(); emp_plus <- c(empresas,"Total")
diag_vol <- data.frame(Serie=character(),Empresa=character(),d_dif=integer(),ADF_p=numeric(),KPSS_p=numeric(),Modelo=character(),RMSE_CV=numeric(),LjungBox_p=numeric(),stringsAsFactors=FALSE)
for (vv in volumenes) {
  th<-rep(0,n_obs); for(emp in empresas){v<-as.numeric((df%>%filter(EMPRESA==emp)%>%arrange(ANNO,TRIMESTRE)%>%pull(!!sym(vv)))[1:n_obs]);v[is.na(v)]<-0;th<-th+v}
  for (emp in emp_plus) {
    key<-paste0(vv,"___",emp); cat("  ->",ETQ[vv],"-",emp,"\n")
    if(emp=="Total"){vals<-th}else{vals<-as.numeric((df%>%filter(EMPRESA==emp)%>%arrange(ANNO,TRIMESTRE)%>%pull(!!sym(vv)))[1:n_obs])}
    vok<-vals[!is.na(vals)&vals!=0]; if(length(vok)<10){res_vol[[key]]<-list(hist=vals,pred=rep(NA,h),lo=rep(NA,h),hi=rep(NA,h),mod="NA");next}
    vals[is.na(vals)]<-0; y_ts<-ts(vals,frequency=4)
    est<-test_estacionariedad(y_ts)
    cv<-tryCatch(cv_temporal(y_ts,h),error=function(e)list(rmse=c(ARIMA_auto=Inf),mejor="ARIMA_auto"))
    tryCatch({
      res<-ajustar_ganador(y_ts,cv$mejor,h); res_vol[[key]]<-list(hist=vals,pred=res$pred,lo=res$lo,hi=res$hi,mod=res$mod)
      diag_vol<-bind_rows(diag_vol,data.frame(Serie=vv,Empresa=emp,d_dif=est$d,ADF_p=round(est$adf_p,4),KPSS_p=round(est$kpss_p,4),Modelo=res$mod,RMSE_CV=round(min(cv$rmse,na.rm=TRUE),2),LjungBox_p=round(res$lb_p,4)))
    },error=function(e){res_vol[[key]]<<-list(hist=vals,pred=rep(tail(vok,1),h),lo=rep(NA,h),hi=rep(NA,h),mod="Naive")})
  }
}

# ═══════════════════════════════════════════════════════════════
# 4. GRUPO C: PRECIOS (competencia + log-transform)
# ═══════════════════════════════════════════════════════════════
cat("[4/7] Grupo C: Precios...\n")
res_precio <- list()
diag_precio <- data.frame(Serie=character(),Empresa=character(),Log=logical(),Modelo=character(),RMSE_CV=numeric(),LjungBox_p=numeric(),stringsAsFactors=FALSE)
for (vp in precios) {
  for (emp in empresas) {
    key<-paste0(vp,"___",emp); cat("  ->",ETQ[vp],"-",emp,"\n")
    vals<-as.numeric((df%>%filter(EMPRESA==emp)%>%arrange(ANNO,TRIMESTRE)%>%pull(!!sym(vp)))[1:n_obs])
    vok<-vals[!is.na(vals)&vals!=0&is.finite(vals)]
    if(length(vok)<10){res_precio[[key]]<-list(hist=vals,pred=rep(NA,h),lo=rep(NA,h),hi=rep(NA,h),mod="NA");next}
    vals[is.na(vals)|is.infinite(vals)]<-0
    usar_log<-all(vok>0)&&(max(vok)/min(vok)>3)
    if(usar_log){y_ts<-ts(log(pmax(vals,1e-6)),frequency=4)}else{y_ts<-ts(vals,frequency=4)}
    cv<-tryCatch(cv_temporal(y_ts,h),error=function(e)list(rmse=c(ETS_damped=Inf),mejor="ETS_damped"))
    tryCatch({
      res<-ajustar_ganador(y_ts,cv$mejor,h)
      if(usar_log){s2<-tryCatch(mean(residuals(ets(y_ts))^2,na.rm=TRUE),error=function(e)0);po<-exp(res$pred+s2/2);lo<-exp(res$lo);hi<-exp(res$hi)
        res_precio[[key]]<-list(hist=vals,pred=po,lo=lo,hi=hi,mod=paste0(res$mod," [log]"))
      }else{res_precio[[key]]<-list(hist=vals,pred=pmax(res$pred,0),lo=pmax(res$lo,0),hi=res$hi,mod=res$mod)}
      diag_precio<-bind_rows(diag_precio,data.frame(Serie=vp,Empresa=emp,Log=usar_log,Modelo=res$mod,RMSE_CV=round(min(cv$rmse,na.rm=TRUE),4),LjungBox_p=round(res$lb_p,4)))
    },error=function(e){res_precio[[key]]<<-list(hist=vals,pred=rep(tail(vok,1),h),lo=rep(NA,h),hi=rep(NA,h),mod="Naive")})
  }
}

# ═══════════════════════════════════════════════════════════════
# 5. GRUPO D: DERIVADOS (HHI, Stenbacka, CR1, Ratio)
# ═══════════════════════════════════════════════════════════════
cat("[5/7] Derivados (HHI, Stenbacka, CR1, Ratio)...\n")

# Funciones identicas al codigo original
calc_hhi <- function(s) 10000 * sum(s^2, na.rm=TRUE)

calc_stenbacka <- function(s) {
  s1 <- max(s, na.rm=TRUE); s2 <- sort(s, decreasing=TRUE)[2]
  g1  <- 0.5 * (1 - (s1^2 - s2^2))
  g05 <- 0.5 * (1 - 0.5 * (s1^2 - s2^2))
  list(sten_g1=g1, sten_g05=g05, s1=s1, dom_g1=ifelse(s1>g1,1,0), dom_g05=ifelse(s1>g05,1,0))
}

derivados <- list()
cp <- c("CUOTA_USUARIOS_VM","CUOTA_TRAFICO_VM","CUOTA_INGRESOS_VM",
        "CUOTA_USUARIOS_IM","CUOTA_TRAFICO_IM","CUOTA_INGRESOS_IM",
        "CUOTA_USUARIOS_SM","CUOTA_INGRESOS_SM")

for (vc in cp) {
  if(is.null(res_cuotas[[vc]])) next; rc <- res_cuotas[[vc]]
  # Historico
  hhi_h <- apply(rc$hist,1,calc_hhi)
  cr1_h <- rc$hist[,"Claro"]*100
  rat_h <- rc$hist[,"Claro"]/pmax(rc$hist[,"Tigo recargado"],1e-6)
  sten_h <- t(sapply(1:n_obs, function(i) { st<-calc_stenbacka(rc$hist[i,]); c(st$sten_g1, st$sten_g05) }))
  # Pronostico
  hhi_p <- apply(rc$pred,1,calc_hhi)
  cr1_p <- rc$pred[,"Claro"]*100
  rat_p <- rc$pred[,"Claro"]/pmax(rc$pred[,"Tigo recargado"],1e-6)
  sten_p <- t(sapply(1:h, function(i) { st<-calc_stenbacka(rc$pred[i,]); c(st$sten_g1, st$sten_g05) }))

  derivados[[vc]] <- list(
    HHI = list(hist=hhi_h, pred=hhi_p),
    CR1 = list(hist=cr1_h, pred=cr1_p),
    Ratio_CR1_CR2 = list(hist=rat_h, pred=rat_p),
    Stenbacka_G1 = list(hist=sten_h[,1]*100, pred=sten_p[,1]*100),
    Stenbacka_G05 = list(hist=sten_h[,2]*100, pred=sten_p[,2]*100)
  )
}

# ═══════════════════════════════════════════════════════════════
# 6. GRAFICAS
# ═══════════════════════════════════════════════════════════════
cat("[6/7] Graficas...\n")

TH <- theme_minimal(base_size=11)+theme(
  plot.background=element_rect(fill="white",color=NA),
  panel.background=element_rect(fill="white",color=NA),
  panel.grid.major=element_line(color="#E0E0E0",linewidth=0.3),panel.grid.minor=element_blank(),
  plot.title=element_text(face="bold",size=14,color="#1E1E1E"),
  plot.subtitle=element_text(size=10,color="#666666"),
  plot.caption=element_text(size=8,color="#A0CD51",hjust=1),
  legend.position="bottom",legend.background=element_rect(fill="white",color=NA),
  axis.text.x=element_text(angle=45,hjust=1,size=7))

grafica_multi <- function(datos_emp, titulo, subtitulo, ylab, fname, es_pct=FALSE) {
  rows <- list()
  for (emp in names(datos_emp)) {
    de <- datos_emp[[emp]]
    if(is.null(de$hist)||all(is.na(de$hist)))next; if(all(is.na(de$pred)))next; if(all(de$hist==0,na.rm=TRUE))next
    ultimo_hist<-de$hist[n_obs]; fecha_corte<-fechas[n_obs]
    rows[[length(rows)+1]]<-data.frame(Fecha=c(fechas,fecha_corte,ff),Valor=c(de$hist,ultimo_hist,de$pred),
      Lower=c(rep(NA,n_obs),NA,de$lo),Upper=c(rep(NA,n_obs),NA,de$hi),
      Tipo=c(rep("H",n_obs),"P",rep("P",h)),Empresa=emp,stringsAsFactors=FALSE)
  }
  if(length(rows)==0)return(invisible(NULL)); dfall<-bind_rows(rows)
  emps_con<-unique(dfall$Empresa); pal_uso<-PAL[emps_con]
  dfall$Empresa<-factor(dfall$Empresa,levels=names(PAL)[names(PAL)%in%emps_con])
  fecha_corte<-fechas[n_obs]; dfall<-dfall%>%filter(!is.na(Valor)&Fecha>=2024)
  if(nrow(dfall)==0)return(invisible(NULL))
  df_labels<-dfall%>%filter(Tipo=="P"&Fecha>fecha_corte)
  df_labels$vj<-ifelse(df_labels$Empresa=="Tigo recargado",2.0,-1.2)
  gg<-ggplot(dfall,aes(x=Fecha,color=Empresa))
  df_ci<-dfall%>%filter(!is.na(Lower)&!is.na(Upper))
  if(nrow(df_ci)>0) gg<-gg+geom_ribbon(data=df_ci,aes(ymin=Lower,ymax=Upper,fill=Empresa),alpha=0.12,color=NA)
  gg<-gg+geom_line(data=dfall%>%filter(Tipo=="H"),aes(y=Valor),linewidth=1.1)+
    geom_line(data=dfall%>%filter(Tipo=="P"),aes(y=Valor),linewidth=1.1,linetype="dashed")+
    geom_point(data=df_labels,aes(y=Valor),size=2.2)+
    geom_vline(xintercept=fecha_corte,linetype="dotted",color="#999999",alpha=0.6)
  if(es_pct){gg<-gg+geom_text(data=df_labels,aes(y=Valor,label=paste0(round(Valor,1),"%"),vjust=vj),size=2.8,show.legend=FALSE)
  }else{gg<-gg+geom_text(data=df_labels,aes(y=Valor,label=scales::comma(round(Valor,0)),vjust=vj),size=2.5,show.legend=FALSE)}
  tf<-sort(unique(dfall$Fecha));tl_all<-c(periodos,pf);tf_all<-c(fechas,ff);tl<-tl_all[match(tf,tf_all)]
  gg<-gg+scale_x_continuous(breaks=tf,labels=tl)+scale_color_manual(values=pal_uso)+scale_fill_manual(values=pal_uso)
  if(es_pct){gg<-gg+scale_y_continuous(labels=function(x)paste0(x,"%"))}else{gg<-gg+scale_y_continuous(labels=scales::comma)}
  gg<-gg+labs(title=titulo,subtitle=subtitulo,x=NULL,y=ylab,color=NULL,fill=NULL,caption="pradatos - Analitica + IA")+TH
  ggsave(file.path(ruta_salida_graficas,fname),gg,width=14,height=7,dpi=200,bg="white")
}

# A: Cuotas
cat("  Graficas A...\n")
for(vc in cuotas){if(is.null(res_cuotas[[vc]]))next;rc<-res_cuotas[[vc]];de<-list()
  for(j in 1:4)de[[empresas[j]]]<-list(hist=rc$hist[,j]*100,pred=rc$pred[,j]*100,lo=rc$lo[,j]*100,hi=rc$hi[,j]*100)
  grafica_multi(de,ETQ[vc],paste0("Modelo: ",rc$mod),"Participacion (%)",paste0("A_",vc,".png"),es_pct=TRUE)}

# B: Volumenes
cat("  Graficas B...\n")
for(vv in volumenes){de<-list();for(emp in emp_plus){key<-paste0(vv,"___",emp);if(is.null(res_vol[[key]]))next;rv<-res_vol[[key]];de[[emp]]<-list(hist=rv$hist,pred=rv$pred,lo=rv$lo,hi=rv$hi)}
  grafica_multi(de,ETQ[vv],"Mejor modelo por CV temporal + IC 95%",NULL,paste0("B_",vv,".png"))}

# C: Precios
cat("  Graficas C...\n")
for(vp in precios){de<-list();for(emp in empresas){key<-paste0(vp,"___",emp);if(is.null(res_precio[[key]]))next;rp<-res_precio[[key]];if(all(is.na(rp$pred)))next
  hv<-rp$hist;hok<-hv[!is.na(hv)&hv!=0&is.finite(hv)];if(length(hok)<5)next;de[[emp]]<-list(hist=rp$hist,pred=rp$pred,lo=rp$lo,hi=rp$hi)}
  grafica_multi(de,ETQ[vp],"Mejor modelo por CV temporal + IC 95%",NULL,paste0("C_",vp,".png"))}

# D: Derivados
cat("  Graficas D...\n")
for(vc in names(derivados)){dd<-derivados[[vc]]
  for(ind in names(dd)){di<-dd[[ind]]
    uh<-di$hist[n_obs];fc0<-fechas[n_obs]
    dfl<-data.frame(Fecha=c(fechas,fc0,ff),Valor=c(di$hist,uh,di$pred),Tipo=c(rep("H",n_obs),"P",rep("P",h)))
    dfl<-dfl%>%filter(Fecha>=2024);df_p<-dfl%>%filter(Tipo=="P"&Fecha>fc0)
    tf<-sort(unique(dfl$Fecha));tl_all<-c(periodos,pf);tf_all<-c(fechas,ff);tl<-tl_all[match(tf,tf_all)]
    gg<-ggplot(dfl,aes(x=Fecha))+
      geom_line(data=dfl%>%filter(Tipo=="H"),aes(y=Valor),color="#A0CD51",linewidth=1.2)+
      geom_line(data=dfl%>%filter(Tipo=="P"),aes(y=Valor),color="#A0CD51",linewidth=1.2,linetype="dashed")+
      geom_point(data=df_p,aes(y=Valor),color="#A0CD51",size=3)+
      geom_text(data=df_p,aes(y=Valor,label=round(Valor,2)),color="#A0CD51",vjust=-1.3,size=3)+
      geom_vline(xintercept=fc0,linetype="dotted",alpha=0.4)+scale_x_continuous(breaks=tf,labels=tl)
    if(ind=="HHI") gg<-gg+geom_hline(yintercept=2500,linetype="dashed",color="#E30613",alpha=0.5)+annotate("text",x=2024,y=2500,label="HHI=2500",size=3,color="#E30613",vjust=-0.5)
    if(ind=="Ratio_CR1_CR2") gg<-gg+geom_hline(yintercept=1.5,linetype="dashed",color="#E30613",alpha=0.5)+annotate("text",x=2024,y=1.5,label="Ratio=1.5",size=3,color="#E30613",vjust=-0.5)
    if(ind=="CR1") gg<-gg+geom_hline(yintercept=50,linetype="dashed",color="#E30613",alpha=0.5)+annotate("text",x=2024,y=50,label="CR1=50%",size=3,color="#E30613",vjust=-0.5)
    gg<-gg+labs(title=paste0(ind," - ",ETQ[vc]),subtitle="Derivado de cuotas pronosticadas",x=NULL,y=ind,caption="pradatos - Analitica + IA")+TH
    ggsave(file.path(ruta_salida_graficas,paste0("D_",ind,"_",vc,".png")),gg,width=12,height=6,dpi=200,bg="white")
  }
}

# ═══════════════════════════════════════════════════════════════
# 7. EXCEL
# ═══════════════════════════════════════════════════════════════
cat("[7/7] Excel...\n")
resumen<-data.frame(Grupo=character(),Serie=character(),Empresa=character(),Modelo=character(),stringsAsFactors=FALSE)
for(vc in names(res_cuotas))resumen<-bind_rows(resumen,data.frame(Grupo="A",Serie=vc,Empresa="Conjunto",Modelo=res_cuotas[[vc]]$mod))
for(key in names(res_vol)){pp<-strsplit(key,"___")[[1]];resumen<-bind_rows(resumen,data.frame(Grupo="B",Serie=pp[1],Empresa=pp[2],Modelo=res_vol[[key]]$mod))}
for(key in names(res_precio)){pp<-strsplit(key,"___")[[1]];resumen<-bind_rows(resumen,data.frame(Grupo="C",Serie=pp[1],Empresa=pp[2],Modelo=res_precio[[key]]$mod))}

hc<-data.frame(Periodo=pf);for(vc in names(res_cuotas)){rc<-res_cuotas[[vc]];for(emp in empresas){cn<-paste0(vc,"_",gsub(" ","_",emp));hc[[cn]]<-round(rc$pred[,emp],6);hc[[paste0(cn,"_LO")]]<-round(rc$lo[,emp],6);hc[[paste0(cn,"_HI")]]<-round(rc$hi[,emp],6)}}
hv<-data.frame(Periodo=pf);for(key in names(res_vol)){rv<-res_vol[[key]];hv[[key]]<-round(rv$pred,2);hv[[paste0(key,"_LO")]]<-round(rv$lo,2);hv[[paste0(key,"_HI")]]<-round(rv$hi,2)}
hp<-data.frame(Periodo=pf);for(key in names(res_precio)){rp<-res_precio[[key]];hp[[key]]<-round(rp$pred,4);hp[[paste0(key,"_LO")]]<-round(rp$lo,4);hp[[paste0(key,"_HI")]]<-round(rp$hi,4)}
hd<-data.frame(Periodo=pf);for(vc in names(derivados))for(ind in names(derivados[[vc]]))hd[[paste0(ind,"_",vc)]]<-round(derivados[[vc]][[ind]]$pred,4)

archivo<-file.path(ruta_salida_excel,"pronosticos_servicios_moviles.xlsx")
writexl::write_xlsx(list(Resumen=resumen,Diag_Volumenes=diag_vol,Diag_Precios=diag_precio,
  A_Cuotas=hc,B_Volumenes=hv,C_Precios=hp,D_Derivados=hd),path=archivo)

cat("\n== COMPLETADO ==\n")
cat("Excel:",archivo,"\n")
cat("Graficas:",ruta_salida_graficas,"\n")
cat("Total:",length(list.files(ruta_salida_graficas,pattern="\\.png$")),"graficas\n")

