# Funciones

## Función para crear flextable
ftable <- function(x) {
  x %>% 
    flextable() %>% 
    theme_vanilla() %>%
    color(part = "footer", color = "#666666") %>%
    color( part = "header", color = "#FFFFFF") %>%
    bg( part = "header", bg = "#2c7fb8") %>%
    fontsize(size = 11) %>%
    font(fontname = 'Calibri') %>%
    # Ajustes de ancho y tipo de alineación de las columnas
    set_table_properties(layout = "autofit") %>% 
    # width(j=1, width = 3) %>%
    align(i = NULL, j = c(2:ncol(x)), align = "right", part = "all")
}

# Cargue y limpieza de datos.

dolar <- read.csv2("datasets/Tasa_de_Cambio_Representativa_del__Mercado_-Historico.csv") %>% 
  clean_names()
hass <- read.xlsx("datasets/Hass_Precios_Historicos.xlsx") %>% clean_names()

# Vamos a sacar un valor promedio de cada variable por mes
#Para el dolar vamos a tomar el campo vigenciadesde

dolar <- dolar %>% 
  mutate(fecha = as.Date(vigenciadesde, format = "%d/%m/%y")) %>% 
  mutate(mes = month(fecha), anio = year(fecha)) %>% 
  mutate(valor = as.double(str_replace(valor,",",""))) %>% 
  group_by(anio,mes) %>%
  summarise(precio_dolar = mean(valor), .groups = "drop") 

# 🖇️ Se genera el vector de meses para usarlo más abajo con la función match.
meses <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo",
           "Junio", "Julio", "Agosto", "Septiembre", "Octubre",
           "Noviembre", "Diciembre")

hass <- hass %>% 
  # Otra forma de sacar el anio.
  #mutate(anio = substring(fecha,nchar(fecha)-4,nchar(fecha))) %>% 
  #mutate(espacio = grepl(", ",fecha))
  mutate(mes = sapply(strsplit(fecha, " "), "[", 2),
         anio = sapply(strsplit(fecha, " "), "[", 5)) %>% 
  mutate(anio = as.double(anio)) %>% 
  # 🖇 Se utiliza la función match con el vector meses.
  mutate(mes = match(mes,meses)) %>% 
  group_by(anio,mes) %>%
  summarise(precio_aguacate_kg = mean(precio_kg), .groups = "drop") 

hass 

hass_dolar <- dolar %>% 
  right_join(hass, by = c("mes","anio")) 

hass_dolar %>% 
  head(5) %>% 
  ftable()

hass_dolar <- hass_dolar %>% mutate(variacion_aguacate = 100*(precio_aguacate_kg -lag(precio_aguacate_kg))/lag(precio_aguacate_kg)) %>% mutate(variacion_dolar = 100*(precio_dolar -lag(precio_dolar))/lag(precio_dolar)) %>% filter(!is.na(variacion_aguacate))

hass_dolar

ggplot(data = hass_dolar,
       aes(x= variacion_dolar, y= variacion_aguacate)) +
geom_point()+ 
geom_hline(yintercept = mean(hass_dolar$variacion_aguacate),
           color = "red")+
geom_text(data=hass_dolar,
          aes(x= 1, y = mean(hass_dolar$variacion_aguacate),
              label = paste("Promedio, [%]:",
                            round(mean(hass_dolar$variacion_aguacate), 2))),
          hjust = 0.5, vjust = -0.1, color = "red"
) +
theme_light()

# 💡 Varianza de los residuales respecto del modelo ingenuo, es la misma varianza de la variable a predecir.

hass_dolar %>% 
  ggplot(aes(x= variacion_dolar, y= variacion_aguacate)) +
  geom_point()+ 
  geom_hline(yintercept = mean(hass_dolar$variacion_aguacate))+
  geom_segment(aes(xend=variacion_dolar, yend=mean(variacion_aguacate)),
               col='red', lty='dashed')+
  theme_light()

#Modelo de regresión lineal (Analítico)

x = hass_dolar$variacion_dolar
y = hass_dolar$variacion_aguacate

# b1 = cov(x,y) / (sigma_x)²
b1 <- cov(x,y) / var(x)
print(b1)

# bo = y_barra - b1*x_barra
b0 <- mean(y) - b1*mean(x)
print(b0)

residuales <- y - (b0 + b1*x)
n <- length(x)

sigma_2 <- sum(residuales^2)/(n-2)
print(sigma_2)

# Confirmación mediante la función lm

mod1 <- hass_dolar %>% lm(variacion_aguacate ~ variacion_dolar,.)
print(mod1$coefficients)

resumen <-  summary(mod1)
(resumen$sigma)^2

parametros <- data.frame(metodo = c("Analítico","Función (lm)"), 
                         Beta_0 = c(b0,mod1$coefficients[1]),
                         Beta_1 = c(b1,mod1$coefficients[2]),
                         sigma_2 = c(sigma_2, (resumen$sigma)^2))

parametros %>% ftable()

tabla_normalidad <- data.frame(
  "Prueba" = c("Kolmogorov-Smirnov (KS)", "Prueba de Lilliefors", "Prueba de Shapiro-Wilk (SW)"),
  "Ventajas" = c(
    "- Puede utilizarse para cualquier distribución teórica.\n- Es una prueba no paramétrica versátil.\n - Adecuada para muestras grandes.",
    "- Es específica para evaluar la normalidad.\n- Puede ser más adecuada para muestras pequeñas.\n - Utiliza estadísticas de prueba que se ajustan a la distribución estándar.",
    "- Es especialmente potente para muestras pequeñas.\n- Diseñada específicamente para evaluar la normalidad.\n- Sensible a desviaciones de la normalidad en cualquier parte de la distribución."
  ),  "Desventajas" = c(
    "- Puede ser menos poderosa en muestras pequeñas.\n - No es específica para la distribución normal.",
    "- Restringida a la comprobación de normalidad. \n - Menos potente en muestras grandes.\n- No es tan versátil como KS para otras distribuciones.",
    "- Más compleja de entender y aplicar.\n- No proporciona estimaciones de parámetros.\n- Puede ser menos adecuada para muestras muy grandes.\n- No es tan versátil como KS para otras distribuciones."))

# Imprimir la tabla con formato

tabla_normalidad %>% 
  ftable() %>% 
  align(i = NULL, j = c(2:3),
        align = "left", part = "all")

# Paso 1: Ordenamos los datos de menor a mayor
n <- length(mod1$residuals)
df_residuals <- data.frame(residuals = mod1$residuals)
df_residuals <- df_residuals %>% 
                arrange(residuals)

# Paso 2: Estimamos los parámetros para la distribución normal teórica a partir de los datos

mean_residuals <- mean(df_residuals$residuals)
sd_residuals <- sd(df_residuals$residuals)

# Paso 3: Calculamos la función acumulada para cada valor de residual
cdf_teorica <- pnorm(df_residuals$residuals,
                      mean = mean_residuals,
                      sd = sd_residuals)

# Paso 4: Calculamos la máxima diferencia absoluta entre la densidad empírica y la teórica

# Método 1: Construir la función empírica usando fórmulas de R (ecdf) y luego tomar diferencias
#           absolutas entre valor empírico y teórico (enfoque KS)
cdf_empirica <- ecdf(df_residuals$residuals)
diferencias_absolutas <- abs(cdf_empirica(df_residuals$residuals) - cdf_teorica)
pos_estadistico_m1 <- which.max(diferencias_absolutas)
estadistico_m1 <- diferencias_absolutas[pos_estadistico_m1]

print(paste0("Este es el valor del estadístico: ", estadistico_m1," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m1))

# Método 2: Construir manualmente los cálculos que permiten hallar la máxima diferencia
#           entre la empírica y la teórica usando la definición de la empírica como
#           función a trozos 📈

e_plus <- seq(1:n)/n
e_minus <- seq(0:(n-1))/n

D_plus <- e_plus - cdf_teorica
D_minus <- cdf_teorica - e_minus

pos_max_d_plus <- which.max(D_plus)
pos_max_d_minus <- which.max(D_minus)

estadistico_m2 <- max(D_plus[pos_max_d_plus], D_minus[pos_max_d_minus])
pos_estadistico_m2 <- function() {
  if(max(D_plus)>=max(D_minus)) {
    return(pos_max_d_plus)
    }else{
      return(pos_max_d_minus)
    }
  }

print(paste0("Este es el valor del estadístico: ", estadistico_m2," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m2()))

# Comprobamos los resultados invocando la función de R que realiza estos cálculos automáticos

Lilliefors_test <- lillie.test(df_residuals$residuals)
print(Lilliefors_test)

estadistico_l <- Lilliefors_test$statistic
print(paste0("Este es el estadístico de Lilliefors usando paquetería de R: ", estadistico_l))

#Gráfico de Lilliefors

x_seg <- df_residuals$residuals[pos_estadistico_m1]
y_sup_seg <- cdf_teorica[pos_estadistico_m1]
y_inf_seg <- cdf_empirica(x_seg)

cat("x_seg: ", format(x_seg, digits=5), 
    " y_inf_seg: ", format(y_inf_seg, digits=5),
    " y_sup_seg: ", format(y_sup_seg, digits=5),
    sep='')

datos_para_grafico <- data.frame(
  residuales = df_residuals$residuals
  )

ggplot(datos_para_grafico, aes(x = residuales)) +
#  geom_density(color = "blue", fill = "lightblue", alpha = 0.5) +
  stat_function(fun = pnorm, args = list(mean = mean_residuals,
                                         sd = sd_residuals),
                color = "red", size = 1) +
  stat_ecdf(color = "green", size = 1) +
  geom_segment(aes(x = x_seg, y = y_inf_seg, xend = x_seg,
                   yend = y_sup_seg),
               color = "purple", size = 1, linetype = "solid") +
  geom_point(aes(x = x_seg, y = y_inf_seg),
             color = "purple",size = 3) +
  geom_point(aes(x = x_seg, y = y_sup_seg),
             color = "purple", size = 3) +
  annotate("text", x = x_seg , y = (y_inf_seg + y_sup_seg)/2, 
           label = paste0("Máxima Separación=",
                          format(estadistico_m1,
                                 digits=5)," 🤔💭", sep=''),
           color = "purple", size = 3, hjust = -0.05, vjust = 0) +
  coord_cartesian(xlim = c(min(datos_para_grafico$residuales),
                           max(datos_para_grafico$residuales))) +
  xlab("Residuales del modelo") +
  ylab("Probabilidad acumulada") +
  labs(title = "Esquema de funcionamiento de la prueba de Lilliefors") +
  theme_minimal()

# Probando normalidad del error. 

shapiro.test(mod1$residuals)

ks.test(mod1$residuals, "pnorm")

# library(lubridate)

hass_dolar$predicciones <- mod1$fitted.values
hass_dolar$residuales <-mod1$residuals
hass_dolar 

# Gráfico de residuales vs valores ajustados

hass_dolar %>% 
  ggplot(aes(x= predicciones, y=residuales)) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales") +
  xlab("Valor predicho")

# Prueba de independencia del error: Durbin-Watson 

resultado_dw <- dwtest(mod1)
print(resultado_dw)

# Prueba de independencia del error: Ljung-Box

p_values <- vector()
for (i in 1:12) {
  resultado_ljung_box <- Box.test(mod1$residuals, lag = i,
                                  type = "Ljung-Box")
  p_values[i] <-resultado_ljung_box$p.value
}

resultados_lb <- data.frame(
  lag = seq(1,12),
  p_value = p_values
)

print(resultado_ljung_box)

print(resultados_lb)

# Residuales estandarizados vs valores predichos

hass_dolar %>% 
  ggplot(aes(x= predicciones, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  xlab("Valor predicho")

# Residuales estandarizados vs precio del dólar

hass_dolar %>% 
  ggplot(aes(x= variacion_dolar, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  xlab("variacion_dolar")

# Pruebas estadísticas formales
# Breusch-Pagan
resultado_breusch_pagan <- bptest(mod1)
print(resultado_breusch_pagan)

# White
resultado_white <- bptest(mod1,
                          varformula = ~ (variacion_dolar +
                                            I(variacion_dolar^2)),
                          data=hass_dolar)
print(resultado_white)

# Implementación manual del test de White

#1. Construir el dataframe que relaciona el cuadrado del error con
#los términos lineales de la predictora y el término cuadrático
#(por ser modelo lineal simple no hay productos cruzados)

df_squared_error <- data.frame(
  squared_error = mod1$residuals ^ 2,
  variacion_dolar = hass_dolar$variacion_dolar,
  variacion_dolar_2 = hass_dolar$variacion_dolar^2
)

#2. Calcule un modelo lineal ajustado sobre el cuadrado del error regresado por el termino lineal y cuadrático del precio del dolar

modelo_white <- lm(data = df_squared_error,
                   squared_error ~ variacion_dolar + variacion_dolar_2)

#3. Calcule el estadístico de prueba como n*R^2 (con n=datos del dataframe de errores)
R2 <- summary(modelo_white)$`r.squared`
w_estadistico = nrow(df_squared_error)*R2
#4. Como el estadístico se distribuye ji-cuadrado con 2 grados de libertad (2 variables regresoras) se tiene que el p valor es:
p_valor = 1-pchisq(as.numeric(w_estadistico),2)
paste("BP = ",round(w_estadistico,4)," y p_valor=",round(p_valor,5))

estimated_squared_error <- modelo_white$coefficients[1] + modelo_white$coefficients[2]*df_squared_error$variacion_dolar +modelo_white$coefficients[3]*df_squared_error$variacion_dolar_2

p<- ggplot(data = df_squared_error, aes(x=variacion_dolar,
                                        y=squared_error))+
  geom_point()+
  geom_line(aes(y=estimated_squared_error, color="red"))+
  labs(title="Modelo white sobre cuadrados del error",
         subtitle = paste("R²=",round(R2,4),
                          "W-Estadístico=",round(w_estadistico,4),
                          "p_valor=",round(p_valor,5)))+
  scale_color_manual(values=c("red"="blue"),
                     labels=c("Error cuadrado estimado"),
                     name="Modelo de White")

print(p)