library(ARDL)
library(forecast)

# 1. Aislar la estacionalidad (Ciclo anual GNSS)
ts_referencia_train <- ts(train_data$height, frequency = 52)
terminos_fourier_train <- fourier(ts_referencia_train, K = 1) 
train_data$sin_anual <- terminos_fourier_train[, 1]
train_data$cos_anual <- terminos_fourier_train[, 2]

# 1. Crear los términos de Fourier
ts_referencia_train <- ts(train_data$height, frequency = 52)
terminos_fourier_train <- fourier(ts_referencia_train, K = 1) 
# 2. Ajustar el modelo estacional para capturar la "Onda" a escala
modelo_estacional <- lm(train_data$height ~ terminos_fourier_train[, 1] + terminos_fourier_train[, 2])
train_data$onda_fourier <- fitted(modelo_estacional)
if (!"fecha" %in% colnames(train_data)) {
  train_data$fecha <- seq(as.Date("2018-07-25"), by = "7 days", length.out = nrow(train_data))
}
ggplot(train_data, aes(x = fecha)) + 
  # Datos originales de fondo (tenues)
  geom_line(aes(y = height, color = "Deformación Bruta"), linewidth = 0.6, alpha = 0.4) +
  # La onda perfecta al frente (gruesa y llamativa)
  geom_line(aes(y = onda_fourier, color = "Ciclo Anual (Fourier K=1)"), linewidth = 1.2) +
  scale_color_manual(values = c("Deformación Bruta" = "#2c3e50", "Ciclo Anual (Fourier K=1)" = "#e74c3c")) +
  theme_bw() +
  labs(
    title = "Aislamiento de la Estacionalidad (Filtro Físico)",
    subtitle = "Extracción de la carga hidrológica anual mediante Serie de Fourier",
    y = "Coordenada Vertical (m)", 
    x = "Fecha de Entrenamiento",
    color = ""
  ) +
  theme(
    legend.position = "bottom", 
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(face = "italic", color = "gray30"),
    legend.text = element_text(face = "bold", size = 11)
  )

# Búsqueda algorítmica del mejor ARDL (Criterio BIC)
seleccion_ardl <- auto_ardl(height ~ east + north | sin_anual + cos_anual, 
                            data = train_data, 
                            max_order = c(4, 4, 4), 
                            selection = "BIC")
mejor_ardl <- seleccion_ardl$best_model
summary(mejor_ardl)

library(sandwich)
library(lmtest)
# 1. Prueba de Cointegración (Bounds Test de Pesaran)
prueba_cointegracion <- bounds_f_test(mejor_ardl, case = 3)
print(prueba_cointegracion)

# 2. Modelo de Corrección de Errores (ECM) con Errores Robustos (HAC)
modelo_ecm <- uecm(mejor_ardl)
ecm_robusto <- coeftest(modelo_ecm, vcov. = NeweyWest(modelo_ecm))

# 3. Multiplicadores de Largo Plazo (Impacto real de la placa tectónica)
multiplicadores_robustos <- multipliers(mejor_ardl, vcov_matrix = NeweyWest(mejor_ardl))
print(multiplicadores_robustos)

# 1. Extraer métricas de precisión
valores_ajustados <- as.numeric(fitted(mejor_ardl))
filas_restantes <- nrow(train_data) - length(valores_ajustados) + 1
valores_reales_train <- train_data$height[filas_restantes:nrow(train_data)]

r_cuadrado <- summary(mejor_ardl)$adj.r.squared

# 2. Gráfica de Ajuste
datos_train_grafica <- train_data[filas_restantes:nrow(train_data), ]
datos_train_grafica$fecha <- fecha_train[filas_restantes:length(fecha_train)]
datos_train_grafica$Pronostico <- valores_ajustados

ggplot(datos_train_grafica, aes(x = fecha)) +
  geom_line(aes(y = height, color = "Altura Real"), linewidth = 0.8, alpha = 0.6) +
  geom_line(aes(y = Pronostico, color = "Ajuste ARDL"), linewidth = 0.8, linetype = "dashed") +
  scale_color_manual(values = c("Altura Real" = "gray60", "Ajuste ARDL" = "darkblue")) +
  theme_minimal() +
  labs(title = "Memoria del Modelo: Ajuste en Datos de Entrenamiento (85%)",
       subtitle = sprintf("Varianza explicada (R-cuadrado Ajustado): %.2f%%", r_cuadrado * 100),
       y = "Desplazamiento Vertical (m)", x = "Fecha", color = "") +
  theme(legend.position = "bottom", plot.title = element_text(face="bold"))

residuos_ardl <- residuals(mejor_ardl)
Box.test(residuos_ardl, type = "Ljung-Box", lag = 10)

lmtest::bgtest(mejor_ardl)

library(dplyr)
library(ggplot2)
library(forecast)

# 1. Proyectar la estacionalidad al futuro
terminos_fourier_test <- fourier(ts_referencia_train, K = 1, h = nrow(test_data))
test_data$sin_anual <- as.numeric(terminos_fourier_test[, 1])
test_data$cos_anual <- as.numeric(terminos_fourier_test[, 2])

# 2. EL ARREGLO: Usar bind_rows en lugar de rbind
data_predict <- bind_rows(train_data, test_data) %>%
  mutate(
    height_L1 = dplyr::lag(height, 1),
    height_L2 = dplyr::lag(height, 2),
    east_L1   = dplyr::lag(east, 1)
  )
train_predict <- data_predict[1:nrow(train_data), ]
test_predict  <- data_predict[(nrow(train_data) + 1):nrow(data_predict), ]
# Re-entrenamos la versión LM explícita con los coeficientes correctos
modelo_lm_predict <- lm(height ~ height_L1 + height_L2 + east + east_L1 + north + sin_anual + cos_anual, 
                        data = train_predict)

# 3. Pronóstico y Cálculo de Error
predicciones_test <- predict(modelo_lm_predict, newdata = test_predict)
rmse_test <- sqrt(mean((test_data$height - predicciones_test)^2, na.rm = TRUE))
mae_test <- mean(abs(test_data$height - predicciones_test), na.rm = TRUE)

# 4. Gráfica Final
# --- SOLUCIÓN: Inyectar fechas al test_data ---
if (!"fecha" %in% colnames(test_data)) {
  # Tomamos la última fecha de tu entrenamiento y le sumamos 7 días
  ultima_fecha_train <- max(train_data$fecha)
  test_data$fecha <- seq(ultima_fecha_train + 7, by = "7 days", length.out = nrow(test_data))
}
df_test_grafica <- data.frame(
  fecha = test_data$fecha, 
  Real = test_data$height,
  Pronostico = as.numeric(predicciones_test)
)

ggplot(df_test_grafica, aes(x = fecha)) +
  geom_line(aes(y = Real, color = "Altura Real"), linewidth = 1) +
  geom_line(aes(y = Pronostico, color = "Pronóstico ARDL"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(values = c("Altura Real" = "gray40", "Pronóstico ARDL" = "darkred")) +
  theme_minimal() +
  labs(title = "Predicción en Datos Nuevos (15%)",
       subtitle = sprintf("Precisión de pronóstico: RMSE = %.5f m | MAE = %.5f m", rmse_test, mae_test),
       y = "Desplazamiento Vertical (m)", x = "Fecha", color = "") +
  theme(legend.position = "bottom", plot.title = element_text(face="bold"))

# RESPUESTA AL IMPULSO (MULTIPLICADORES DINÁMICOS) PARA ARDL
coefs <- coef(modelo_lm_predict)
horizonte <- 20
respuestas <- numeric(horizonte)

# Simulamos un "Shock" unitario en el componente ESTE
# Recordando la estructura: Height ~ Height_L1 + Height_L2 + East + East_L1 ...
impulse_value <- 1 

# Calculamos el efecto inmediato y propagado (Recursión)
for (t in 1:horizonte) {
  if (t == 1) {
    # Efecto inmediato: Coeficiente de 'east'
    respuestas[t] <- coefs["east"] * impulse_value
  } else if (t == 2) {
    # Primer rezago: Coeficiente de 'east_L1' + persistencia de la altura
    respuestas[t] <- (coefs["height_L1"] * respuestas[t-1]) + (coefs["east_L1"] * impulse_value)
  } else {
    # Efecto a largo plazo: Solo queda la persistencia (memoria) de la altura
    respuestas[t] <- (coefs["height_L1"] * respuestas[t-1]) + (coefs["height_L2"] * respuestas[t-2])
  }
}

# 4. Graficar el Multiplicador Dinámico
df_irf_ardl <- data.frame(Semana = 1:horizonte, Impacto = respuestas)

ggplot(df_irf_ardl, aes(x = Semana, y = Impacto)) +
  geom_line(color = "blue", linewidth = 1) +
  geom_point(color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Respuesta Dinámica de la Altura ante un Shock en el Este (ARDL)",
       subtitle = "Muestra cómo el hundimiento se estabiliza tras el movimiento horizontal",
       x = "Semanas después del evento", y = "Impacto acumulado en Altura (m)")

# 1. Preparar la matriz XREG para el futuro (Test Data)
# IMPORTANTE: Usamos los datos de prueba que el modelo NUNCA ha visto
xreg_test_arimax <- cbind(
  east = test_data$east,
  north = test_data$north,
  sin_anual = test_data$sin_anual,
  cos_anual = test_data$cos_anual
)
pronostico_arimax <- forecast(modelo_sarimax, xreg = xreg_test_arimax)

rmse_arimax <- sqrt(mean((test_data$height - pronostico_arimax$mean)^2, na.rm = TRUE))
mae_arimax <- mean(abs(test_data$height - pronostico_arimax$mean), na.rm = TRUE)

df_grafica_arimax <- data.frame(
  fecha = test_data$fecha,
  Real = test_data$height,
  Pronostico = as.numeric(pronostico_arimax$mean)
)
checkresiduals(modelo_sarimax)

ggplot(df_grafica_arimax, aes(x = fecha)) +
  geom_line(aes(y = Real, color = "Altura Real"), linewidth = 1) +
  geom_line(aes(y = Pronostico, color = "Pronóstico ARIMAX"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(values = c("Altura Real" = "gray40", "Pronóstico ARIMAX" = "blue")) +
  theme_minimal() +
  labs(title = "Validación Externa: Pronóstico ARIMAX (Out-of-Sample)",
       subtitle = sprintf("Precisión: RMSE = %.5f m | MAE = %.5f m", rmse_arimax, mae_arimax),
       y = "Desplazamiento Vertical (m)", x = "Fecha", color = "") +
  theme(legend.position = "bottom", plot.title = element_text(face="bold"))

library(vars)
library(urca)
matriz_gnss <- cbind(
  height = train_data$height,
  east = train_data$east,
  north = train_data$north
)
seleccion_lags <- VARselect(matriz_gnss, lag.max = 10, type = "const")
#cat("\n--- REZAGOS ÓPTIMOS SUGERIDOS ---\n")
#print(seleccion_lags$selection)
K_optimo <- max(seleccion_lags$selection["SC(n)"], 2)
cat("\nRezago (K) elegido para Johansen:", K_optimo, "\n")

# type = "trace": Usamos el estadístico de Traza (el más robusto)
# ecdet = "const": Asumimos intercepto pero sin tendencia en la ecuación (equivalente a tu Case 3)
prueba_johansen <- ca.jo(matriz_gnss, type = "trace", ecdet = "const", K = K_optimo)
cat("\n--- RESULTADOS DE LA PRUEBA DE JOHANSEN ---\n")

summary(prueba_johansen)

# ==============================================================================
#  CONSTRUCCIÓN DEL MODELO Y SIMULACIÓN DE CHOQUES
# ==============================================================================
# 1. Convertimos el test de Johansen en un modelo calculable
# Usamos r = 1 (asumiendo 1 relación principal de equilibrio a pesar del over-rejection)
modelo_vecm <- vec2var(prueba_johansen, r = 1)

# 2. Funciones de Impulso-Respuesta (IRF)
# Simulamos un "shock" en el Este y vemos cómo reacciona la Altura 20 periodos adelante
irf_este_altura <- irf(modelo_vecm, impulse = "east", response = "height", 
                       n.ahead = 20, boot = TRUE)

# Al correr esto, R te abrirá una ventana con una gráfica física espectacular
plot(irf_este_altura, main = "Impacto en la Altura tras un 'Shock' hacia el Este")

# ==============================================================================
# FASE 5 (VECM): PREDICCIÓN DEL FUTURO Y CÁLCULO DEL RMSE
# ==============================================================================
# 3. Le pedimos al modelo que pronostique tantos pasos adelante como tenga tu test_data
prediccion_vecm <- predict(modelo_vecm, n.ahead = nrow(test_data))

# Extraemos exclusivamente los valores pronosticados para la Altura
pronostico_altura_vecm <- prediccion_vecm$fcst$height[, "fcst"]

# Calculamos el error final (RMSE y MAE)
rmse_vecm <- sqrt(mean((test_data$height - pronostico_altura_vecm)^2, na.rm = TRUE))
mae_vecm <- mean(abs(test_data$height - pronostico_altura_vecm), na.rm = TRUE)
cat("RMSE VECM:", round(rmse_vecm, 5), "metros\n")

cat("MAE VECM:", round(mae_vecm, 5), "metros\n")

v_real   <- as.numeric(test_data$height)
v_ardl   <- as.numeric(predicciones_test)
v_arimax <- as.numeric(pronostico_arimax$mean)
v_vecm   <- as.numeric(pronostico_altura_vecm)
df_integrado <- data.frame(
  fecha  = test_data$fecha, 
  Real   = v_real,
  ARDL   = v_ardl,
  ARIMAX = v_arimax,
  VECM   = v_vecm,
  ARDL_lwr = v_ardl - (1.96 * rmse_test),
  ARDL_upr = v_ardl + (1.96 * rmse_test)
)
ggplot(df_integrado, aes(x = fecha)) +
  geom_ribbon(aes(ymin = ARDL_lwr, ymax = ARDL_upr, fill = "Incertidumbre ARDL (95%)"), alpha = 0.15) +
  geom_line(aes(y = Real, color = "Observado (Real)"), linewidth = 1.1) +
  geom_line(aes(y = ARDL, color = "ARDL (Dinámico)"), linewidth = 0.9) +
  geom_line(aes(y = ARIMAX, color = "ARIMAX-Fourier"), linewidth = 0.8, linetype = "dashed") +
  geom_line(aes(y = VECM, color = "VECM (Sistémico)"), linewidth = 0.8, linetype = "dotdash") +
  scale_color_manual(values = c(
    "Observado (Real)" = "black", 
    "ARDL (Dinámico)"  = "blue", 
    "ARIMAX-Fourier"   = "darkgreen", 
    "VECM (Sistémico)" = "darkred"
  )) +
  scale_fill_manual(values = c("Incertidumbre ARDL (95%)" = "dodgerblue")) +
  theme_minimal() +
  labs(title = "Comparación de Modelos: Predicción de Altura Wanchaq",
       subtitle = "Evaluación sobre el 15% de datos finales (no vistos)",
       y = "Desplazamiento Vertical (m)", x = "Año",
       color = "Modelos", fill = "") +
  theme(legend.position = "bottom", plot.title = element_text(face="bold", size = 14))

# Función para el Nash-Sutcliffe Efficiency
f_nse <- function(real, pred) {
  1 - (sum((real - pred)^2, na.rm=T) / sum((real - mean(real, na.rm=T))^2, na.rm=T))
}

# Consolidación de la tabla
tabla_final <- data.frame(
  Modelo = c("ARDL", "ARIMAX-Fourier", "VECM"),
  RMSE_m = c(rmse_test, rmse_arimax, rmse_vecm), 
  MAE_m  = c(mae_test, mae_arimax, mae_vecm),
  NSE    = c(f_nse(test_data$height, predicciones_test),
             f_nse(test_data$height, as.numeric(pronostico_arimax$mean)),
             f_nse(test_data$height, as.numeric(pronostico_altura_vecm)))
)
# Redondeo para estética académica
tabla_final[, 2:3] <- round(tabla_final[, 2:3], 5)
tabla_final[, 4]   <- round(tabla_final[, 4], 4)

print(tabla_final)