# 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")
}

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

dolar %>% 
  head(10)%>% 
  ftable()

hass%>% 
  head(10)%>% 
  ftable() 

# 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_dolar <- dolar %>% 
  right_join(hass, by = c("mes","anio")) 

hass_dolar %>% 
  head(10) %>% 
  ftable()

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


# dolar<- dolar%>% 
#   mutate(variacion_dolar= (precio_dolar- lag(precio_dolar))/lag(precio_dolar)) %>% 
#   filter(!is.na(variacion_dolar))

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

hass_dolar %>% 
  ggplot(aes(x=variacion_dolar, y= variacion_aguacate )) +
  geom_point()+ theme_light()

media_aguacate <- mean(hass_dolar$variacion_aguacate)

ggplot(data = hass_dolar, aes(x = variacion_dolar, y = variacion_aguacate)) +
  geom_point() +
  geom_hline(yintercept = media_aguacate, color = "red", linetype = "dashed") +
  annotate("text", x = max(hass_dolar$variacion_dolar), y = media_aguacate, 
           label = paste("Valor medio del Aguacate:", round(media_aguacate, 4)),
           hjust = 1, vjust = 0, color = "red") +
  labs(title = "Modelo Ingenuo", x = "Variación del Dólar", y = "Variación del Aguacate") +
  theme_minimal()

n <- length(hass_dolar$variacion_aguacate)
MSE_ingenuo  <- (1/n) * sum((hass_dolar$variacion_aguacate - mean(hass_dolar$variacion_aguacate))^2)
MSE_ingenuo

hass_dolar %>%
  ggplot(aes(x = variacion_dolar, y = variacion_aguacate)) +
  geom_point() +
  geom_hline(yintercept = mean(hass_dolar$variacion_aguacate), color = "red") +
  geom_segment(aes(xend = variacion_dolar, yend = mean(hass_dolar$variacion_aguacate)), col = 'red', linetype = 'dashed') +
  labs(title = "Distancia entre cada dato observado y el promedio histórico", x = "Capital", y = "FNCER") +
  theme_light()

b0 <- 0.002
b1 <- 0.007



# Predicciones
x <- hass_dolar$variacion_dolar
y_pred <- b0 + b1 * x

# Gráfica

hass_dolar %>%
  ggplot(aes(x =variacion_dolar , y = variacion_aguacate)) +
  geom_point() +
  geom_line(aes(x = x, y = y_pred), col = "red") +
  labs(x = "Variación_Dolar", y = "Variación_aguacate") +
  theme_light()

# Extraemos las variables de interés
x = hass_dolar$variacion_dolar
y = hass_dolar$variacion_aguacate

# Calculamos medias muestrales para las variables X y Y

x_barra <- mean(x)
y_barra <- mean(y)

# Definimos posibles valores para b1, y por consiguiente para b0
b1_values <- seq(-1, 0.5, by = 0.08)
b0_values <- y_barra - b1_values * x_barra

# Fijamos límites en X y Y para cada gráfico
xmin = min(x) - 0
xmax = max(x) + 1
ymin = min(y) - 0
ymax = max(y) + 1

# Construimos una función para calcular el MSE en función de x_i, y_i, b1 y b0
MSE <- function(x, y, b0, b1) {
  errors <- y - (b0 + b1*x)
  squared_errors <- errors^2
  mse <- mean(squared_errors)
  return(mse)
}

# Crear un dataframe para almacenar los resultados
list_mse <- vector()

# Calcular el MSE para diferentes valores de B1
for (i in 1:length(b1_values)) {
  mse <- MSE(x, y, b0_values[i], b1_values[i])
  list_mse[i] <- mse
}

mse_df = data.frame(
  B0 = b0_values,
  B1 = b1_values,
  mse = list_mse,
  color_point = 'green'
)

# Capturando el renglón donde ocurre el mínimo mse
pos_min <- which(mse_df$mse==min(mse_df$mse))

# Impriendo las filas cercanas al mínimo
mse_df$color_point[(pos_min - 4):(pos_min + 4)] <- rep('red',9)
ftable(mse_df[(pos_min - 8):(pos_min + 8), ])

if (!file.exists("beta_vs_mse.gif")){
  # Crear la animación
  animation <- ggplot(mse_df, aes(B1, mse)) +
    geom_line(aes( x= B1, y = mse, color = color_point), linetype = "dashed") +
    geom_point(aes(color = color_point), size = 3) +
    labs(title = "Valores de MSE para cada posible B1, respetando E(e) = 0") +
    geom_text(data = mse_df, 
              aes(label = paste("MSE: ", round(mse,0))), x= 0.5, y = 0.5, color='black'
    ) +
    geom_text(data = mse_df, 
              aes(label = paste("MSE mínimo: ", round(min(mse),0))),
              x= 1, y = min(mse_df$mse) - 1e5, color='red') +
    scale_color_manual(values = c("red" = "red", "green" = "green"),
                       labels = c("Subóptima", "Óptima"),
                       name = "Soluciones") +
    theme_minimal() +
    transition_reveal(B1)

  animate(animation, duration = 10, fps = 2, renderer = gifski_renderer())
  anim_save("beta_vs_mse.gif")
}

if (!file.exists("beta_estimation.gif")){
  
# Las pausas en la animación son fotogramas con información repetida correspondiente
# al renglón óptimo

fotogramas_pausa <- data.frame(
  B0 = rep(b0_values[pos_min], 2*nrow(mse_df)),
  B1 = rep(b1_values[pos_min], 2*nrow(mse_df)),
  mse = rep(list_mse[pos_min], 2*nrow(mse_df)),
  color_point = rep('red',2*nrow(mse_df))  
)

mse_df <- rbind(mse_df, fotogramas_pausa)
mse_df <- mse_df[order(mse_df$B1), ]
mse_df$estado <- seq(1,nrow(mse_df))

# Construimos el dataframe de predicciones para cada estado, con el cual poder
# graficar los residuales en cada fotograma

df_predicciones <- data.frame(B0_pred=numeric(0),
                              B1_pred=numeric(0),
                              x_from_pred=numeric(0),
                              y_from_pred=numeric(0),
                              y_pred_mod=numeric(0),
                              estado = numeric(0))
for (i in 1:nrow(mse_df)){
  df_nuevo <-data.frame(
    BO_pred = rep(mse_df[i, 'B0'], length(x)),
    B1_pred = rep(mse_df[i, 'B1'], length(x)),
    x_from_pred = x,
    y_from_pred = y,
    y_pred_mod = mse_df[i, 'B0'] + mse_df[i,'B1']*x,
    estado = rep(mse_df[i, 'estado'], length(x))
  ) 
  df_predicciones <- rbind(df_predicciones, df_nuevo)
}
w <- data.frame(df_predicciones)
# Creamos el gráfico base

base_plot <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
  geom_point() +
  geom_point(x=mean(x), y=mean(y), color='green', size=2.5) +
  geom_abline(aes(slope = 0, intercept = y_barra), 
              linetype = "dashed", color = "blue") +
  labs(title = "Estimación de Betas por Minimización del MSE") +
  theme_minimal()

# Creamos la animación

animation <- base_plot +
  geom_abline(aes(slope = B1, intercept = B0, color = color_point), data=mse_df) +
  geom_text(aes(label = paste("MSE: ", round(mse,0))), x= 0.5, y = 0.5, data=mse_df,
            color = "red")+
  geom_segment(data = df_predicciones,
               aes(x=x_from_pred, y=y_from_pred, xend=x_from_pred, yend=y_pred_mod),
               col = "violet", lty='dashed') +
  transition_states(estado, transition_length = 2, state_length = 1) +
  scale_color_manual(values = c("red" = "red", "green" = "green"),
                     labels = c("Subóptima", "Óptima"),
                     name = "Rectas regresoras") +
  enter_fade() +
  exit_fade()

# Guardamos la animación en un archivo GIF
animate(animation, duration = 20, fps = 2, renderer = gifski_renderer())
anim_save("beta_estimation.gif", animation)
}

# 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)

# Confirmamos usando la función lm, a continuación la documentación

# lm(formula, data, subset, weights, na.action,
#    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
#    singular.ok = TRUE, contrasts = NULL, offset, ...)

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

resumen <-  summary(mod1)
resumen$sigma

parametros <- data.frame(metodo = c("manual","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()

# 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))

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")

# Primero inspeccionamos las variables individuales para observar su
# comportamiento temporal en un mismo gráfico

# library(lubridate)
x <- hass_dolar$variacion_dolar
y <- hass_dolar$variacion_aguacate


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

# Gráfico de residuales vs valores ajustados

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

# Gráfico de residuales con índice temporal
hass_dolar %>% 
  ggplot(aes(x= variacion_dolar, y=variacion_aguacate)) +
  geom_point()+ 
  geom_line(aes(y = residuales, color = "residuales"), size = 1)+
  theme_light()+
  xlab("capital")+
  ylab("Residuales")

# Prueba de Durbin-Watson tomando el mismo orden en el que venían los
# datos en el dataframe puesto que este orden se corresponde con el 
# orden temporal.

resultado_dw <- dwtest(mod1)
print(resultado_dw)

# Prueba de Ljung-Box
# Como esperamos cierta estacionalidad intra año en los datos podría tener sentido esperarla también para los residuales. Por lo tanto eligiremos lag = 1,2,3,...,12

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 fncer")

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

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

resultado_white <- bptest(mod1,
  varformula = ~ (hass_dolar$variacion_dolar + I(hass_dolar$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)

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = hass_dolar$variacion_dolar,
  capital_2 = hass_dolar$variacion_dolar^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el estadístico de prueba como n * R^2 (con n = número de datos del dataframe de errores)
R2 <- summary(modelo_white)$r.squared
n <- nrow(df_squared_error)
w_estadistico <- n * R2

# 4. Calcular el valor p utilizando la distribución chi-cuadrado con 2 grados de libertad
p_valor <- 1 - pchisq(w_estadistico, df = 2)

# Imprimir el estadístico de prueba y el valor p
paste("BP =", round(w_estadistico, 4), "y p_valor =", round(p_valor, 5))

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = hass_dolar$variacion_dolar,
  capital_2 = hass_dolar$variacion_dolar^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el valor estimado del error cuadrado
estimated_squared_error <- modelo_white$coefficients[1] +
                          modelo_white$coefficients[2] * df_squared_error$capital +
                          modelo_white$coefficients[3] * df_squared_error$capital_2

# 4. Crear un gráfico
library(ggplot2)

p <- ggplot(data = df_squared_error, aes(x = capital, 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)

# 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)

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = hass_dolar$variacion_dolar,
  capital_2 = hass_dolar$variacion_dolar^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el estadístico de prueba como n * R^2 (con n = número de datos del dataframe de errores)
R2 <- summary(modelo_white)$r.squared
n <- nrow(df_squared_error)
w_estadistico <- n * R2

# 4. Calcular el valor p utilizando la distribución chi-cuadrado con 2 grados de libertad
p_valor <- 1 - pchisq(w_estadistico, df = 2)

# Imprimir el estadístico de prueba y el valor p
paste("BP =", round(w_estadistico, 4), "y p_valor =", round(p_valor, 5))

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = hass_dolar$variacion_dolar,
  capital_2 = hass_dolar$variacion_dolar^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el valor estimado del error cuadrado
estimated_squared_error <- modelo_white$coefficients[1] +
                          modelo_white$coefficients[2] * df_squared_error$capital +
                          modelo_white$coefficients[3] * df_squared_error$capital_2

# 4. Crear un gráfico
library(ggplot2)

p <- ggplot(data = df_squared_error, aes(x = capital, 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)