# Extración de Residuos
residuos <- residuals(modelo_seguro)
options(scipen = 99) # Evitar notación científica en los resultados
# A) Prueba de Shapiro-Wilk
# Calculo de W
salida_SW <- shapiro.test(modelo_seguro$residuals)
print(salida_SW)

# Resultado final de Wn
Wn_salida<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(Wn_salida)

# B) Prueba de Jarque-Bera
library(tseries)
options(scipen = 99) # Evitar notación científica en los resultados
jb_test <- jarque.bera.test(residuos)
print(jb_test)

options(scipen = 99) # Evitar notación científica en los resultados
# C) Prueba de Kolmogorov-Smirnov
library(nortest)
prueba_KS<-lillie.test(modelo_seguro$residuals)
prueba_KS

library(fastGraph) 
# Parámetros del nivel de significancia y estadísticos
alpha_sig <- 0.05
JB_stat <- jb_test$statistic
gl <- jb_test$parameter

# Cálculo del Valor Crítico Teórico
VC_jb <- qchisq(1 - alpha_sig, gl, lower.tail = TRUE)

# Grafico
shadeDist(VC_jb, 
          ddist = "dchisq", 
          parm1 = gl, 
          lower.tail = FALSE, 
          xmin = 0, 
          xmax = 10, 
          sub = paste("VC:", round(VC_jb, 2), " ", "JB:", round(JB_stat, 2)))

# OPCION 1 
#  Cálculo del Indice de Condición usando librería “mctest”
library(mctest)
X_mat<-model.matrix(modelo_seguro)
mctest(mod = modelo_seguro)

# CÁLCULO DE LAS VIFs usando performance
library(performance)
VIFs<-multicollinearity(x = modelo_seguro,verbose = FALSE)
VIFs

# CÁLCULO DE LAS VIFs usando mctest
library(mctest)
mc.plot(mod = modelo_seguro,vif = 2)

library(lmtest)
options(scipen = 99) # Evitar notación científica en los resultados
# Ejecutar la prueba de White para el nuevo modelo
PWhite <- bptest(
  modelo_seguro,
  ~ I(age^2) + I(bmi^2) + I(children^2) + age * bmi + age * children + age * smoker + age * region + bmi * children + bmi * smoker + bmi * region +  children * smoker + children * region + smoker * region, data = insurance)

# Mostrar los resultados
print(PWhite)

## Usando Libreraria Skedastic (Otra alternativa)
library(skedastic)
options(scipen = 99) # Evitar notación científica en los resultados
skedastic::white(modelo_seguro, interactions = FALSE) 

library(lmtest)
dwtest(modelo_seguro,alternative = "two.sided",iterations = 1000)

library(car)
durbinWatsonTest(modelo_seguro,simulate = TRUE,reps = 1000)

library(lmtest)
bgtest(modelo_seguro,order = 2)

library(lmtest)
bgtest(modelo_seguro,order = 1)

# Extraemos los residuos de tu modelo
residuos <- residuals(modelo_seguro)

# Graficamos con ZOOM
acf(residuos, 
    lag.max = 5,                      # Zoom en X: Muestra solo hasta el lag 5
    ylim = c(-0.2, 0.2),              # Zoom en Y: Corta el gráfico entre -0.2 y 0.2
    main = "Correlograma de los Residuos",
    xlab = "Lag", 
    ylab = "ACF")

library(stargazer)
options(scipen = 2)
#Data para la predicción X'm basado en un perfil de cliente para insurance
X_m<-data.frame(age=35, bmi=28.5, children=2, smoker="no", region="southwest")
# Intervalos de Confianza del 95% y del 99%
confidense<-c(0.95,0.99)
#Predicción usando predict (calculando ambos niveles para la matriz)
pred_95 <- predict(object = modelo_seguro, newdata = X_m, interval = "prediction", level = 0.95)
pred_99 <- predict(object = modelo_seguro, newdata = X_m, interval = "prediction", level = 0.99)

# Reconvertir de log a dólares aplicando exp()
pred_95_exp <- exp(pred_95)
pred_99_exp <- exp(pred_99)

predicciones <- list()
predicciones$fit <- rbind(pred_95_exp, pred_99_exp)

rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
          title = "Pronósticos e intervalos de confianza (en dólares)",
          type = "html") #Poner results='asis' en opciones del chunk

library(forecast)
library(kableExtra)
options(scipen = 2)
#Data para la predicción X'm
X_m<-data.frame(age=35, bmi=28.5, children=2, smoker="no", region="southwest")
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95,0.99)

#Realizando el pronóstico con forecast
pronosticos<-forecast(object = modelo_seguro,
         level = confidense,
         newdata = X_m, ts = FALSE)

# Reconvertir de log a dólares aplicando exp()
pronosticos_df <- as.data.frame(pronosticos)
pronosticos_exp <- exp(pronosticos_df)

kable(pronosticos_exp,
      caption = "Pronóstico e intervalos de confianza (en dólares):",
      digits = 1, format = "html") #Poner results='asis' en opciones del chunk

#Bias Proportion
Um<-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado) 
}
#Variance Proportion
Us<-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc<-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U<-function(pronosticado,observado){
   library(DescTools)
  RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}

options(scipen = 99) #No mostrar notación cientifica.
library(dplyr) # Para manejo de datos y activar el operador "pipe" %>%
library(caret) # Permite Realizar muestreo sobre los data frame
library(DescTools) # Contiene las funciones para calcular las medidas de performance
library(stargazer) # Para dar formato, y obtener resumen estadistico de las simulaciones
set.seed(50) # Permite fijar la semilla aleatoria, para reproducir los resultados obtenidos en esta clase

# SOLICITADO: Simulación con 5,000 réplicas / iteraciones
numero_de_muestras<-5000 

# Se crea la lista con las 5000 muestras sobre los cargos de seguros (80% entrenamiento)
muestras<- insurance$charges %>%
  createDataPartition(p = 0.8,
                      times = numero_de_muestras,
                      list = TRUE)

# Listas vacias, que contendran los datos de entrenamiento, los pronosticos para los datos de prueba, y para las estadisticas de cada muestra
Modelos_Entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance<-vector(mode = "list",
                              length = numero_de_muestras)

#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- insurance[muestras[[j]], ]
Datos_Prueba<- insurance[-muestras[[j]], ]

# Modelo adaptado a tu estructura justificada
Modelos_Entrenamiento[[j]]<-lm(formula = log(charges) ~ age + bmi + children + smoker + region, data = Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)

Resultados_Performance_data_entrenamiento[[j]]<-data.frame( 
            R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
                    log(Datos_Entrenamiento$charges)),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        log(Datos_Entrenamiento$charges)),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      log(Datos_Entrenamiento$charges)),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       log(Datos_Entrenamiento$charges))*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         log(Datos_Entrenamiento$charges),type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges)),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges)),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges))
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            RMSE = RMSE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            MAE = MAE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            MAPE= MAPE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges))*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], log(Datos_Prueba$charges),
                         type = 1), # También se puede usar la función que creamos: THEIL_U
            Um=Um(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            Us=Us(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            Uc=Uc(Pronostico_Prueba[[j]], log(Datos_Prueba$charges))
            )
} #No olvidar este corchete ;)

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo (Entrenamiento 80%)",
            type = "html",
            digits = 3)

bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación (Predictivo - 5,000 Réplicas)",
            type = "html",
            digits = 3)