Trabajo final econometría 1

#SECCIÓN 0: IMPORTAR LAS LIBRERIAS

# install.packages(c("readxl", "dplyr", "tidyr", "ggplot2", "lmtest",
#                    "tseries", "sandwich", "car", "nortest", "stargazer",
#                    "scales", "gridExtra"))

library(readxl)    # Leer archivos Excel (.xlsx)

## Warning: package 'readxl' was built under R version 4.4.3

library(dplyr)     # Manipular datos: filtrar, calcular columnas, unir tablas

## Warning: package 'dplyr' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)     # Reorganizar datos (formato ancho → largo, útil para gráficas)

## Warning: package 'tidyr' was built under R version 4.4.3

library(ggplot2)   # Crear gráficas profesionales

## Warning: package 'ggplot2' was built under R version 4.4.3

library(lmtest)    # Pruebas diagnósticas: Breusch-Pagan, Durbin-Watson, etc.

## Warning: package 'lmtest' was built under R version 4.4.3

## Cargando paquete requerido: zoo

## 
## Adjuntando el paquete: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tseries)   # Prueba de raíz unitaria Dickey-Fuller Aumentado (adf.test)

## Warning: package 'tseries' was built under R version 4.4.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(sandwich)  # Errores estándar robustos tipo Newey-West

## Warning: package 'sandwich' was built under R version 4.4.3

library(car)       # VIF: detectar multicolinealidad entre regresores

## Warning: package 'car' was built under R version 4.4.3

## Cargando paquete requerido: carData

## 
## Adjuntando el paquete: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(nortest)   # Prueba de normalidad Lilliefors
library(stargazer) # Exportar tablas de resultados en formato académico

## 
## Please cite as:

##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

library(scales)    # Formatear ejes (porcentajes, separadores de miles, etc.)

## Warning: package 'scales' was built under R version 4.4.3

library(gridExtra) # Combinar múltiples gráficas en una sola figura

## Warning: package 'gridExtra' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

SECCIÓN 1: RUTAS Y CARGA DE DATOS

# Ruta del archivo Excel con la base consolidada
ruta_base <- "C:/Users/diego/OneDrive/Imagenes1/Escritorio/univalle/8 SEMESTRE/ECONOMETRIA/TRABAJO FINAL/Base_Final_Econometria.xlsx"

# Carpeta donde se guardarán todas las gráficas y tablas generadas
ruta_salida <- "C:/Users/diego/OneDrive/Imagenes1/Escritorio/univalle/8 SEMESTRE/ECONOMETRIA/TRABAJO FINAL/"

cat("Cargando base de datos...\n")

## Cargando base de datos...

datos_raw <- read_excel(ruta_base, sheet = "Sheet1")

# Renombramos las columnas con nombres claros y consistentes
colnames(datos_raw) <- c("fecha_num", "anio", "trimestre",
                         "desempleo", "subempleo", "ocupacion", "PET",
                         "hurto_per", "hurto_com", "homicidios",
                         "nacimientos", "defunciones", "IPC")

# Creamos una fecha real en formato Date para las gráficas.
# El mes de inicio de cada trimestre: Q1=enero(1), Q2=abril(4),
# Q3=julio(7), Q4=octubre(10)
datos <- datos_raw %>%
  mutate(
    mes_ini = case_when(
      trimestre == 1 ~ 1,
      trimestre == 2 ~ 4,
      trimestre == 3 ~ 7,
      trimestre == 4 ~ 10
    ),
    fecha = as.Date(paste(anio, mes_ini, 1, sep = "-")),
    t     = row_number()   # Índice temporal secuencial: 1, 2, 3, ..., 39
  ) %>%
  arrange(fecha)

cat(paste("  Observaciones:", nrow(datos), "\n"))

##   Observaciones: 39

cat(paste("  Variables:", ncol(datos), "\n"))

##   Variables: 16

cat(paste("  Período:", format(min(datos$fecha)), "a", format(max(datos$fecha)), "\n"))

##   Período: 2015-01-01 a 2024-07-01

cat(paste("  Valores faltantes (NA):", sum(is.na(datos[, 4:13])), "\n\n"))

##   Valores faltantes (NA): 0

SECCIÓN 2: ESTADÍSTICAS DESCRIPTIVAS

# Las estadísticas descriptivas son el primer paso del análisis empírico.
# Permiten entender el comportamiento de cada variable antes de modelar:
#   - Media: valor promedio en todo el período
#   - Desv. Estándar: qué tanto varían los datos alrededor de la media
#   - Mínimo / Máximo: valores extremos (útiles para detectar outliers)
#   - Mediana: valor central (menos sensible a outliers que la media)
# ==============================================================================

cat("=== ESTADÍSTICAS DESCRIPTIVAS ===\n\n")

## === ESTADÍSTICAS DESCRIPTIVAS ===

vars_desc <- datos %>%
  select(desempleo, subempleo, ocupacion, PET,
         hurto_per, hurto_com, homicidios,
         nacimientos, defunciones, IPC)

# Tabla descriptiva manual con las estadísticas más relevantes
tabla_desc <- data.frame(
  Variable  = c("Tasa Desempleo", "Tasa Subempleo", "Tasa Ocupación",
                "Pob. Edad Trabajar (miles)", "Hurto a Personas",
                "Hurto al Comercio", "Homicidios",
                "Nacimientos", "Defunciones", "IPC (var. mensual prom.)"),
  Media     = round(sapply(vars_desc, mean,   na.rm = TRUE), 4),
  Desv.Est  = round(sapply(vars_desc, sd,     na.rm = TRUE), 4),
  Minimo    = round(sapply(vars_desc, min,    na.rm = TRUE), 4),
  Mediana   = round(sapply(vars_desc, median, na.rm = TRUE), 4),
  Maximo    = round(sapply(vars_desc, max,    na.rm = TRUE), 4),
  row.names = NULL
)

print(tabla_desc)

##                      Variable     Media Desv.Est    Minimo   Mediana    Maximo
## 1              Tasa Desempleo    0.1337   0.0388    0.0939    0.1246    0.2982
## 2              Tasa Subempleo    0.0976   0.0371    0.0000    0.1083    0.1341
## 3              Tasa Ocupación    0.5644   0.0421    0.4008    0.5760    0.6187
## 4  Pob. Edad Trabajar (miles) 1793.1538  60.0564 1700.0000 1797.0000 1894.0000
## 5            Hurto a Personas 1425.4530 495.3016  607.3333 1423.6667 2634.3333
## 6           Hurto al Comercio  223.6154 103.5672   64.0000  227.6667  454.3333
## 7                  Homicidios   95.8120  14.9861   66.3333   95.3333  134.6667
## 8                 Nacimientos 5739.2051 870.2847 4157.0000 5812.0000 7250.0000
## 9                 Defunciones 3793.0256 652.2821 3235.0000 3524.0000 5928.0000
## 10   IPC (var. mensual prom.)    0.0048   0.0040   -0.0010    0.0042    0.0154

# Exportar en formato texto (para copiar al documento del trabajo)
stargazer(as.data.frame(vars_desc),
          type    = "text",
          title   = "Tabla 1. Estadísticas Descriptivas – Cali 2015 Q1 a 2024 Q3",
          digits  = 4,
          out     = paste0(ruta_salida, "Tabla1_Descriptivas.txt"))

## 
## Tabla 1. Estadísticas Descriptivas – Cali 2015 Q1 a 2024 Q3
## ======================================================
## Statistic   N     Mean    St. Dev.   Min       Max    
## ------------------------------------------------------
## desempleo   39   0.1337    0.0388   0.0939    0.2982  
## subempleo   39   0.0976    0.0371   0.0000    0.1341  
## ocupacion   39   0.5644    0.0421   0.4008    0.6187  
## PET         39 1,793.1540 60.0564   1,700     1,894   
## hurto_per   39 1,425.4530 495.3016 607.3333 2,634.3330
## hurto_com   39  223.6154  103.5672 64.0000   454.3333 
## homicidios  39  95.8120   14.9861  66.3333   134.6667 
## nacimientos 39 5,739.2050 870.2847  4,157     7,250   
## defunciones 39 3,793.0260 652.2821  3,235     5,928   
## IPC         39   0.0048    0.0040  -0.0010    0.0154  
## ------------------------------------------------------

cat("\nTabla 1 exportada como Tabla1_Descriptivas.txt\n\n")

## 
## Tabla 1 exportada como Tabla1_Descriptivas.txt

SECCIÓN 3: GRÁFICAS EN NIVELES

# Antes de modelar, siempre hay que visualizar las series en el tiempo.
# Buscamos identificar:
#   1. TENDENCIA: ¿la serie sube o baja de forma sostenida con el tiempo?
#      Ejemplo: PET sube cada año porque Cali crece → tiene tendencia.
#   2. QUIEBRES ESTRUCTURALES: ¿hay saltos bruscos?
#      El COVID-19 en 2020 generó el pico máximo de desempleo.
#   3. OUTLIERS: valores extremos alejados del patrón general.
#
# Esta inspección visual nos ayuda a decidir qué transformaciones
# son necesarias antes de estimar el modelo.
# ==============================================================================

cat("Generando gráficas de series en niveles...\n")

## Generando gráficas de series en niveles...

# --- Gráfica 1: Variable dependiente – Tasa de Desempleo ---
g1 <- ggplot(datos, aes(x = fecha, y = desempleo)) +
  geom_line(color = "#C0392B", linewidth = 1.1) +
  geom_point(size = 2, color = "#C0392B") +
  # Línea vertical para marcar el quiebre estructural del COVID-19
  geom_vline(xintercept = as.Date("2020-03-01"),
             linetype = "dashed", color = "gray40", linewidth = 0.8) +
  annotate("text", x = as.Date("2020-06-01"),
           y = max(datos$desempleo, na.rm = TRUE) * 0.98,
           label = "COVID-19", size = 3.2, color = "gray30", hjust = 0) +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(
    title    = "Gráfica 1. Tasa de Desempleo en Cali (2015 Q1 – 2024 Q3)",
    subtitle = "Datos trimestrales fijos | Fuente: DANE – GEIH / Cali Cómo Vamos",
    x = NULL,
    y = "Tasa de Desempleo"
  ) +
  theme_minimal(base_size = 12)

ggsave(paste0(ruta_salida, "Grafica1_Desempleo.png"),
       g1, width = 10, height = 5, dpi = 150)
print(g1)

# --- Gráfica 2: Variables de seguridad (los 3 regresores clave) ---
# pivot_longer convierte 3 columnas en formato largo:
# en lugar de 3 columnas separadas, tenemos 1 columna "variable" y 1 "valor"
# Esto permite usar facet_wrap para hacer un panel con 3 gráficas individuales
seg_long <- datos %>%
  select(fecha, homicidios, hurto_per, hurto_com) %>%
  pivot_longer(cols = -fecha, names_to = "variable", values_to = "valor") %>%
  mutate(variable = case_when(
    variable == "homicidios" ~ "Homicidios",
    variable == "hurto_per"  ~ "Hurto a Personas",
    variable == "hurto_com"  ~ "Hurto al Comercio",
    TRUE ~ variable
  ))

g2 <- ggplot(seg_long, aes(x = fecha, y = valor, color = variable)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = as.Date("2020-03-01"),
             linetype = "dashed", color = "gray40", linewidth = 0.7) +
  facet_wrap(~variable, scales = "free_y", ncol = 1) +
  labs(
    title    = "Gráfica 2. Indicadores de Inseguridad en Cali (2015 Q1 – 2024 Q3)",
    subtitle = "Promedios trimestrales | Fuente: Observatorio de Seguridad / Policía Nacional",
    x = NULL,
    y = "Promedio trimestral de eventos"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

ggsave(paste0(ruta_salida, "Grafica2_Seguridad.png"),
       g2, width = 10, height = 9, dpi = 150)
print(g2)

# --- Gráfica 3: Variables del mercado laboral ---
ml_long <- datos %>%
  select(fecha, desempleo, subempleo, ocupacion) %>%
  pivot_longer(cols = -fecha, names_to = "variable", values_to = "valor") %>%
  mutate(variable = case_when(
    variable == "desempleo" ~ "Tasa de Desempleo",
    variable == "subempleo" ~ "Tasa de Subempleo",
    variable == "ocupacion" ~ "Tasa de Ocupación",
    TRUE ~ variable
  ))

g3 <- ggplot(ml_long, aes(x = fecha, y = valor, color = variable)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = as.Date("2020-03-01"),
             linetype = "dashed", color = "gray40", linewidth = 0.7) +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  facet_wrap(~variable, scales = "free_y", ncol = 1) +
  labs(
    title    = "Gráfica 3. Variables del Mercado Laboral en Cali (2015 Q1 – 2024 Q3)",
    subtitle = "Fuente: DANE – GEIH / Cali Cómo Vamos",
    x = NULL,
    y = "Tasa"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

ggsave(paste0(ruta_salida, "Grafica3_MercadoLaboral.png"),
       g3, width = 10, height = 9, dpi = 150)
print(g3)

cat("Gráficas 1, 2 y 3 exportadas.\n\n")

## Gráficas 1, 2 y 3 exportadas.

SECCIÓN 4: PRUEBA DE RAÍZ UNITARIA (Dickey-Fuller Aumentado – ADF)

# ¿QUÉ ES UNA RAÍZ UNITARIA?
# Una serie tiene raíz unitaria cuando NO es estacionaria: su media y/o
# varianza cambian con el tiempo porque tiene una tendencia persistente.
# Ejemplos: PIB, precios, población. En tu base, el PET sube año tras año.
#
# ¿POR QUÉ ES UN PROBLEMA EN REGRESIÓN?
# Si regresas una variable con tendencia creciente sobre otra con tendencia
# creciente, el modelo va a mostrar una relación "significativa" aunque
# en realidad no exista ninguna causalidad real. Eso se llama REGRESIÓN
# ESPURIA: R² alto y coeficientes significativos, pero estadísticamente falso.
#
# CÓMO FUNCIONA LA PRUEBA ADF:
#   H0 (hipótesis nula): la serie TIENE raíz unitaria → NO es estacionaria
#   H1 (alternativa):    la serie NO tiene raíz unitaria → SÍ es estacionaria
#
#   Regla de decisión:
#   p-valor < 0.05 → rechazamos H0 → la serie ES estacionaria ✓
#   p-valor ≥ 0.05 → NO rechazamos H0 → hay raíz unitaria → DIFERENCIAR
# ==============================================================================

cat("=== PRUEBA ADF – SERIES EN NIVELES ===\n")

## === PRUEBA ADF – SERIES EN NIVELES ===

cat("H0: la serie tiene raíz unitaria (no es estacionaria)\n")

## H0: la serie tiene raíz unitaria (no es estacionaria)

cat("Regla: p < 0.05 → estacionaria ✓  |  p ≥ 0.05 → raíz unitaria → diferenciar\n\n")

## Regla: p < 0.05 → estacionaria ✓  |  p ≥ 0.05 → raíz unitaria → diferenciar

vars_nivel <- c("desempleo", "subempleo", "ocupacion", "PET",
                "hurto_per", "hurto_com", "homicidios",
                "nacimientos", "defunciones", "IPC")

nombres_nivel <- c("Desempleo", "Subempleo", "Ocupación", "PET",
                   "Hurto Personas", "Hurto Comercio", "Homicidios",
                   "Nacimientos", "Defunciones", "IPC")

tabla_adf_nivel <- data.frame(
  Variable    = character(),
  Estadistico = numeric(),
  p_valor     = numeric(),
  Conclusion  = character(),
  stringsAsFactors = FALSE
)

for (i in seq_along(vars_nivel)) {
  serie  <- na.omit(datos[[ vars_nivel[i] ]])
  result <- adf.test(serie)
  tabla_adf_nivel <- rbind(tabla_adf_nivel, data.frame(
    Variable    = nombres_nivel[i],
    Estadistico = round(result$statistic, 3),
    p_valor     = round(result$p.value,   4),
    Conclusion  = ifelse(result$p.value < 0.05,
                         "Estacionaria",
                         "Raiz unitaria → Diferenciar"),
    stringsAsFactors = FALSE
  ))
}

print(tabla_adf_nivel, row.names = FALSE)

##        Variable Estadistico p_valor                  Conclusion
##       Desempleo      -1.625  0.7203 Raiz unitaria → Diferenciar
##       Subempleo      -1.970  0.5849 Raiz unitaria → Diferenciar
##       Ocupación      -1.512  0.7646 Raiz unitaria → Diferenciar
##             PET      -2.822  0.2515 Raiz unitaria → Diferenciar
##  Hurto Personas      -2.358  0.4333 Raiz unitaria → Diferenciar
##  Hurto Comercio      -1.290  0.8513 Raiz unitaria → Diferenciar
##      Homicidios      -3.184  0.1098 Raiz unitaria → Diferenciar
##     Nacimientos      -1.635  0.7164 Raiz unitaria → Diferenciar
##     Defunciones      -2.646  0.3203 Raiz unitaria → Diferenciar
##             IPC      -1.474  0.7795 Raiz unitaria → Diferenciar

write.csv(tabla_adf_nivel,
          paste0(ruta_salida, "Tabla2_ADF_Niveles.csv"),
          row.names = FALSE)
cat("\nTabla 2 exportada como Tabla2_ADF_Niveles.csv\n\n")

## 
## Tabla 2 exportada como Tabla2_ADF_Niveles.csv

SECCIÓN 5: TRANSFORMACIÓN – PRIMERAS DIFERENCIAS

# ¿QUÉ ES LA PRIMERA DIFERENCIA?
#   Δy_t = y_t − y_{t−1}
#
# En lugar de trabajar con el NIVEL de cada variable (ej: el desempleo
# fue 12% este trimestre), trabajamos con el CAMBIO entre períodos
# (ej: el desempleo SUBIÓ 1.5 puntos respecto al trimestre anterior).
#
# ¿POR QUÉ ESTO ELIMINA LA TENDENCIA?
# Si una variable sube 0.5 puntos cada trimestre de forma constante,
# la diferencia siempre será ≈ 0.5: un número sin dirección → estacionaria.
#
# COSTO: perdemos 1 observación (la primera fila queda NA porque no
# tiene período anterior), por eso pasamos de 39 a 38 observaciones.
#
# INTERPRETACIÓN DE LOS COEFICIENTES:
# β₁ mide: "cuando los homicidios CAMBIAN en 1 unidad respecto al
# trimestre anterior, ¿cuánto CAMBIA la tasa de desempleo ese trimestre?"
# ==============================================================================

cat("=== CONSTRUYENDO PRIMERAS DIFERENCIAS ===\n\n")

## === CONSTRUYENDO PRIMERAS DIFERENCIAS ===

# lag(x) devuelve el valor del período ANTERIOR de x
# Ejemplo: desempleo   = c(0.12, 0.14, 0.11)
#          lag(desempleo) = c(NA,  0.12, 0.14)
#          Δdesempleo   = c(NA, +0.02, -0.03)

datos_d <- datos %>%
  arrange(fecha) %>%
  mutate(
    d_desempleo   = desempleo   - lag(desempleo),   # Variable DEPENDIENTE
    d_subempleo   = subempleo   - lag(subempleo),
    d_ocupacion   = ocupacion   - lag(ocupacion),
    d_PET         = PET         - lag(PET),
    d_hurto_per   = hurto_per   - lag(hurto_per),
    d_hurto_com   = hurto_com   - lag(hurto_com),
    d_homicidios  = homicidios  - lag(homicidios),
    d_nacimientos = nacimientos - lag(nacimientos),
    d_defunciones = defunciones - lag(defunciones),
    d_IPC         = IPC         - lag(IPC)
  ) %>%
  filter(!is.na(d_desempleo))   # Eliminar la primera fila (NA inevitable)

cat(paste("  Observaciones después de diferenciar:", nrow(datos_d),
          "(39 − 1 = 38)\n\n"))

##   Observaciones después de diferenciar: 38 (39 − 1 = 38)

# Verificar que las series diferenciadas ya son estacionarias
cat("=== PRUEBA ADF – SERIES DIFERENCIADAS ===\n")

## === PRUEBA ADF – SERIES DIFERENCIADAS ===

cat("Esperamos p < 0.05 en todas (confirmación de estacionariedad)\n\n")

## Esperamos p < 0.05 en todas (confirmación de estacionariedad)

vars_d     <- c("d_desempleo", "d_subempleo", "d_ocupacion", "d_PET",
                "d_hurto_per", "d_hurto_com", "d_homicidios",
                "d_nacimientos", "d_defunciones", "d_IPC")

nombres_d  <- c("Δ Desempleo", "Δ Subempleo", "Δ Ocupación", "Δ PET",
                "Δ Hurto Personas", "Δ Hurto Comercio", "Δ Homicidios",
                "Δ Nacimientos", "Δ Defunciones", "Δ IPC")

tabla_adf_dif <- data.frame(
  Variable   = character(),
  p_valor    = numeric(),
  Conclusion = character(),
  stringsAsFactors = FALSE
)

for (i in seq_along(vars_d)) {
  serie  <- na.omit(datos_d[[ vars_d[i] ]])
  result <- adf.test(serie)
  tabla_adf_dif <- rbind(tabla_adf_dif, data.frame(
    Variable   = nombres_d[i],
    p_valor    = round(result$p.value, 4),
    Conclusion = ifelse(result$p.value < 0.05,
                        "Estacionaria ✓",
                        "Aun no estacionaria"),
    stringsAsFactors = FALSE
  ))
}

## Warning in adf.test(serie): p-value smaller than printed p-value
## Warning in adf.test(serie): p-value smaller than printed p-value

print(tabla_adf_dif, row.names = FALSE)

##          Variable p_valor          Conclusion
##       Δ Desempleo  0.2056 Aun no estacionaria
##       Δ Subempleo  0.0100      Estacionaria ✓
##       Δ Ocupación  0.2232 Aun no estacionaria
##             Δ PET  0.0827 Aun no estacionaria
##  Δ Hurto Personas  0.0580 Aun no estacionaria
##  Δ Hurto Comercio  0.0100      Estacionaria ✓
##      Δ Homicidios  0.0693 Aun no estacionaria
##     Δ Nacimientos  0.0102      Estacionaria ✓
##     Δ Defunciones  0.5297 Aun no estacionaria
##             Δ IPC  0.0942 Aun no estacionaria

write.csv(tabla_adf_dif,
          paste0(ruta_salida, "Tabla3_ADF_Diferencias.csv"),
          row.names = FALSE)
cat("\nTabla 3 exportada como Tabla3_ADF_Diferencias.csv\n\n")

## 
## Tabla 3 exportada como Tabla3_ADF_Diferencias.csv

# --- Gráfica 4: Nivel vs Diferencia del desempleo ---
# Esta gráfica es muy útil para el trabajo: muestra visualmente por qué
# hay que diferenciar y qué efecto tiene la transformación
g4a <- ggplot(datos, aes(x = fecha, y = desempleo)) +
  geom_line(color = "#C0392B", linewidth = 1) +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(title = "En NIVELES (con tendencia visible)",
       x = NULL, y = "Nivel") +
  theme_minimal(base_size = 11)

g4b <- ggplot(datos_d, aes(x = fecha, y = d_desempleo)) +
  geom_line(color = "#2980B9", linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_y_continuous(labels = percent_format(accuracy = 0.1)) +
  labs(title = "En PRIMERAS DIFERENCIAS (estacionaria)",
       x = NULL, y = "Cambio trimestral") +
  theme_minimal(base_size = 11)

g4 <- grid.arrange(g4a, g4b, ncol = 1,
                   top = "Gráfica 4. Efecto de Diferenciar la Tasa de Desempleo")

ggsave(paste0(ruta_salida, "Grafica4_Diferenciacion.png"),
       g4, width = 10, height = 7, dpi = 150)
cat("Gráfica 4 exportada.\n\n")

## Gráfica 4 exportada.

SECCIÓN 6: ESTIMACIÓN DEL MODELO OLS EN PRIMERAS DIFERENCIAS

# MODELO ECONOMÉTRICO COMPLETO:
#
# Δdesempleo_t = β₀ + β₁·Δhomicidios_t + β₂·Δhurto_personas_t
#              + β₃·Δhurto_comercio_t  + β₄·Δsubempleo_t
#              + β₅·Δocupacion_t       + β₆·ΔPET_t
#              + β₇·Δdefunciones_t     + β₈·Δnacimientos_t
#              + β₉·ΔIPC_t + ε_t
#
# OLS (Mínimos Cuadrados Ordinarios):
# Encuentra los valores de β que minimizan la suma de cuadrados de los
# residuos (Σε²): es decir, la línea que más se acerca a los datos.
#
# CÓMO LEER EL SUMMARY():
#   Estimate (β̂):  efecto marginal estimado
#   Std. Error:    incertidumbre en la estimación
#   t value:       β̂ / Std.Error (qué tan lejos está β̂ de cero)
#   Pr(>|t|):      p-valor. Si < 0.05 → significativo al 5%
#                           Si < 0.10 → significativo al 10%
#   R² ajustado:   % de la varianza de Δdesempleo explicada por el modelo
#   F-estadístico: prueba si el modelo en conjunto es significativo
# ==============================================================================

cat("=== ESTIMACIÓN OLS – MODELO COMPLETO ===\n\n")

## === ESTIMACIÓN OLS – MODELO COMPLETO ===

modelo1 <- lm(
  d_desempleo ~ d_homicidios + d_hurto_per + d_hurto_com +
    d_subempleo  + d_ocupacion + d_PET        +
    d_defunciones + d_nacimientos + d_IPC,
  data = datos_d
)

print(summary(modelo1))

## 
## Call:
## lm(formula = d_desempleo ~ d_homicidios + d_hurto_per + d_hurto_com + 
##     d_subempleo + d_ocupacion + d_PET + d_defunciones + d_nacimientos + 
##     d_IPC, data = datos_d)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0253280 -0.0061087 -0.0002461  0.0078156  0.0244521 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.250e-03  6.336e-03   0.355  0.72512    
## d_homicidios  -8.064e-06  1.778e-04  -0.045  0.96415    
## d_hurto_per   -1.480e-05  1.087e-05  -1.362  0.18391    
## d_hurto_com   -1.581e-06  4.718e-05  -0.034  0.97350    
## d_subempleo    1.872e-01  6.766e-02   2.767  0.00991 ** 
## d_ocupacion   -9.891e-01  9.472e-02 -10.442 3.67e-11 ***
## d_PET         -4.511e-04  1.175e-03  -0.384  0.70386    
## d_defunciones  8.442e-06  5.286e-06   1.597  0.12153    
## d_nacimientos  2.290e-06  6.610e-06   0.346  0.73165    
## d_IPC          1.072e+00  5.590e-01   1.918  0.06534 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01281 on 28 degrees of freedom
## Multiple R-squared:  0.8838, Adjusted R-squared:  0.8465 
## F-statistic: 23.67 on 9 and 28 DF,  p-value: 8.074e-11

stargazer(modelo1,
          type             = "text",
          title            = "Tabla 4. Modelo OLS Completo en Primeras Diferencias – Cali 2015‑2024",
          dep.var.labels   = "Δ Tasa de Desempleo",
          covariate.labels = c("Δ Homicidios", "Δ Hurto Personas", "Δ Hurto Comercio",
                               "Δ Subempleo",  "Δ Ocupación",      "Δ PET",
                               "Δ Defunciones","Δ Nacimientos",    "Δ IPC"),
          out              = paste0(ruta_salida, "Tabla4_Modelo_Completo.txt"))

## 
## Tabla 4. Modelo OLS Completo en Primeras Diferencias – Cali 2015‑2024
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                         Δ Tasa de Desempleo    
## -----------------------------------------------
## Δ Homicidios                 -0.00001          
##                              (0.0002)          
##                                                
## Δ Hurto Personas             -0.00001          
##                              (0.00001)         
##                                                
## Δ Hurto Comercio             -0.00000          
##                              (0.00005)         
##                                                
## Δ Subempleo                  0.187***          
##                               (0.068)          
##                                                
## Δ Ocupación                  -0.989***         
##                               (0.095)          
##                                                
## Δ PET                         -0.0005          
##                               (0.001)          
##                                                
## Δ Defunciones                 0.00001          
##                              (0.00001)         
##                                                
## Δ Nacimientos                 0.00000          
##                              (0.00001)         
##                                                
## Δ IPC                         1.072*           
##                               (0.559)          
##                                                
## Constant                       0.002           
##                               (0.006)          
##                                                
## -----------------------------------------------
## Observations                    38             
## R2                             0.884           
## Adjusted R2                    0.846           
## Residual Std. Error       0.013 (df = 28)      
## F Statistic           23.667*** (df = 9; 28)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

cat("\nTabla 4 exportada.\n\n")

## 
## Tabla 4 exportada.

# ------------------------------------------------------------------------------
# MODELO PARSIMONIOSO
# Con base en los resultados del modelo completo, construimos un modelo más
# simple que incluye:
#   (a) Las variables de seguridad SIEMPRE (son el núcleo de la hipótesis)
#   (b) Las variables de control que resultaron significativas (p < 0.10)
#
# *** AJUSTA las variables de abajo según los p-valores del modelo1 ***
# Por defecto incluimos: variables de seguridad + ocupación + IPC
# ------------------------------------------------------------------------------

cat("=== MODELO PARSIMONIOSO ===\n")

## === MODELO PARSIMONIOSO ===

cat("(Seguridad + variables significativas del modelo completo)\n\n")

## (Seguridad + variables significativas del modelo completo)

modelo2 <- lm(
  d_desempleo ~ d_homicidios + d_hurto_per + d_hurto_com +
    d_ocupacion  + d_IPC,
  data = datos_d
)

print(summary(modelo2))

## 
## Call:
## lm(formula = d_desempleo ~ d_homicidios + d_hurto_per + d_hurto_com + 
##     d_ocupacion + d_IPC, data = datos_d)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.035319 -0.009270  0.000217  0.008623  0.036759 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.005e-04  2.284e-03  -0.044   0.9652    
## d_homicidios  1.716e-04  1.724e-04   0.995   0.3270    
## d_hurto_per  -1.163e-05  1.128e-05  -1.031   0.3104    
## d_hurto_com  -3.072e-05  4.515e-05  -0.680   0.5012    
## d_ocupacion  -9.214e-01  1.001e-01  -9.205 1.65e-10 ***
## d_IPC         1.376e+00  5.266e-01   2.613   0.0136 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01393 on 32 degrees of freedom
## Multiple R-squared:  0.8429, Adjusted R-squared:  0.8184 
## F-statistic: 34.34 on 5 and 32 DF,  p-value: 5.823e-12

# Tabla comparativa de ambos modelos en una sola hoja
stargazer(modelo1, modelo2,
          type             = "text",
          title            = "Tabla 5. Comparación de Modelos OLS en Primeras Diferencias",
          dep.var.labels   = "Δ Tasa de Desempleo",
          column.labels    = c("Modelo Completo", "Parsimonioso"),
          covariate.labels = c("Δ Homicidios", "Δ Hurto Personas", "Δ Hurto Comercio",
                               "Δ Subempleo",  "Δ Ocupación",      "Δ PET",
                               "Δ Defunciones","Δ Nacimientos",    "Δ IPC"),
          out              = paste0(ruta_salida, "Tabla5_Comparacion_Modelos.txt"))

## 
## Tabla 5. Comparación de Modelos OLS en Primeras Diferencias
## =================================================================
##                                  Dependent variable:             
##                     ---------------------------------------------
##                                  Δ Tasa de Desempleo             
##                        Modelo Completo          Parsimonioso     
##                              (1)                    (2)          
## -----------------------------------------------------------------
## Δ Homicidios               -0.00001                0.0002        
##                            (0.0002)               (0.0002)       
##                                                                  
## Δ Hurto Personas           -0.00001               -0.00001       
##                           (0.00001)              (0.00001)       
##                                                                  
## Δ Hurto Comercio           -0.00000               -0.00003       
##                           (0.00005)              (0.00005)       
##                                                                  
## Δ Subempleo                0.187***                              
##                            (0.068)                               
##                                                                  
## Δ Ocupación               -0.989***              -0.921***       
##                            (0.095)                (0.100)        
##                                                                  
## Δ PET                      -0.0005                               
##                            (0.001)                               
##                                                                  
## Δ Defunciones              0.00001                               
##                           (0.00001)                              
##                                                                  
## Δ Nacimientos              0.00000                               
##                           (0.00001)                              
##                                                                  
## Δ IPC                       1.072*                1.376**        
##                            (0.559)                (0.527)        
##                                                                  
## Constant                    0.002                 -0.0001        
##                            (0.006)                (0.002)        
##                                                                  
## -----------------------------------------------------------------
## Observations                  38                     38          
## R2                          0.884                  0.843         
## Adjusted R2                 0.846                  0.818         
## Residual Std. Error    0.013 (df = 28)        0.014 (df = 32)    
## F Statistic         23.667*** (df = 9; 28) 34.344*** (df = 5; 32)
## =================================================================
## Note:                                 *p<0.1; **p<0.05; ***p<0.01

cat("\nTabla 5 exportada.\n\n")

## 
## Tabla 5 exportada.

SECCIÓN 7: DIAGNÓSTICOS DEL MODELO

# Para que las inferencias (pruebas t, intervalos de confianza) sean válidas,
# los residuos del modelo deben cumplir 4 supuestos de Gauss-Markov.
# Verificamos cada uno con una prueba estadística formal y una gráfica.
# ==============================================================================

cat("=== DIAGNÓSTICOS DEL MODELO COMPLETO ===\n\n")

## === DIAGNÓSTICOS DEL MODELO COMPLETO ===

residuos  <- residuals(modelo1)
ajustados <- fitted(modelo1)
datos_d$residuos  <- residuos
datos_d$ajustados <- ajustados

# ------------------------------------------------------------------
# SUPUESTO 1 – NORMALIDAD DE RESIDUOS
#
# Los errores deben seguir una distribución normal. Sin normalidad,
# las pruebas t y F pueden no ser válidas en muestras pequeñas (n < 50).
#
# Shapiro-Wilk:  H0 = los residuos son normales
# Lilliefors:    alternativa más potente para n moderado
#   p > 0.05 → no rechazamos H0 → hay normalidad ✓
#   p < 0.05 → rechazamos H0   → no hay normalidad (problema)
# ------------------------------------------------------------------

cat("--- Supuesto 1: Normalidad de residuos ---\n")

## --- Supuesto 1: Normalidad de residuos ---

sw <- shapiro.test(residuos)
ll <- lillie.test(residuos)

cat(paste("  Shapiro-Wilk:  W =", round(sw$statistic, 4),
          "| p-valor =", round(sw$p.value, 4), "|",
          ifelse(sw$p.value > 0.05, "Normalidad ✓", "No normalidad"), "\n"))

##   Shapiro-Wilk:  W = 0.9821 | p-valor = 0.7893 | Normalidad ✓

cat(paste("  Lilliefors:    D =", round(ll$statistic, 4),
          "| p-valor =", round(ll$p.value, 4), "|",
          ifelse(ll$p.value > 0.05, "Normalidad ✓", "No normalidad"), "\n\n"))

##   Lilliefors:    D = 0.0877 | p-valor = 0.6515 | Normalidad ✓

# Histograma de residuos con curva normal teórica superpuesta
g5 <- ggplot(data.frame(res = residuos), aes(x = res)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12,
                 fill = "#3498DB", color = "white", alpha = 0.85) +
  stat_function(fun  = dnorm,
                args = list(mean = mean(residuos), sd = sd(residuos)),
                color = "red", linewidth = 1.2) +
  labs(title    = "Gráfica 5. Distribución de los Residuos",
       subtitle = "La curva roja es la distribución normal teórica",
       x = "Residuos", y = "Densidad") +
  theme_minimal(base_size = 12)

ggsave(paste0(ruta_salida, "Grafica5_Histograma_Residuos.png"),
       g5, width = 8, height = 5, dpi = 150)
print(g5)

# Q-Q Plot: si los puntos siguen la línea diagonal → normalidad
png(paste0(ruta_salida, "Grafica6_QQplot.png"),
    width = 800, height = 600, res = 150)
qqnorm(residuos,
       main = "Gráfica 6. Q-Q Plot de Residuos",
       xlab = "Cuantiles teóricos (normal)",
       ylab = "Cuantiles de los residuos")
qqline(residuos, col = "red", lwd = 2)
dev.off()

## png 
##   2

cat("Gráficas 5 y 6 exportadas.\n\n")

## Gráficas 5 y 6 exportadas.

# ------------------------------------------------------------------
# SUPUESTO 2 – HOMOCEDASTICIDAD
#
# La varianza de los errores debe ser CONSTANTE para todas las
# observaciones. Si cambia (HETEROCEDASTICIDAD), los errores
# estándar de los coeficientes son incorrectos → p-valores falsos.
#
# Prueba Breusch-Pagan:
#   H0: varianza constante (homocedasticidad)
#   p < 0.05 → heterocedasticidad → usar errores estándar ROBUSTOS
# ------------------------------------------------------------------

cat("--- Supuesto 2: Homocedasticidad (Breusch-Pagan) ---\n")

## --- Supuesto 2: Homocedasticidad (Breusch-Pagan) ---

bp <- bptest(modelo1)
cat(paste("  BP =", round(bp$statistic, 4),
          "| p-valor =", round(bp$p.value, 4), "|",
          ifelse(bp$p.value > 0.05,
                 "Homocedasticidad ✓",
                 "Heterocedasticidad → usar errores robustos"), "\n\n"))

##   BP = 7.274 | p-valor = 0.6086 | Homocedasticidad ✓

# Errores estándar robustos Newey-West:
# Corrigen simultáneamente heterocedasticidad Y autocorrelación.
# Si hay cualquiera de los dos problemas, estos errores son más confiables.
cat("  Coeficientes con errores estándar robustos (Newey-West, lag = 2):\n")

##   Coeficientes con errores estándar robustos (Newey-West, lag = 2):

nw <- coeftest(modelo1, vcov = NeweyWest(modelo1, lag = 2))
print(round(nw, 5))

## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.00225    0.00258  0.8708  0.39124    
## d_homicidios  -0.00001    0.00008 -0.1003  0.92085    
## d_hurto_per   -0.00001    0.00001 -1.5623  0.12946    
## d_hurto_com    0.00000    0.00003 -0.0491  0.96122    
## d_subempleo    0.18721    0.03839  4.8760    4e-05 ***
## d_ocupacion   -0.98911    0.10633 -9.3026  < 2e-16 ***
## d_PET         -0.00045    0.00046 -0.9746  0.33813    
## d_defunciones  0.00001    0.00000  2.4331  0.02161 *  
## d_nacimientos  0.00000    0.00000  0.4862  0.63060    
## d_IPC          1.07236    0.40914  2.6210  0.01401 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n")

# Gráfica residuos vs valores ajustados:
# Debe verse como una nube aleatoria alrededor de 0.
# Un patrón en forma de abanico o embudo indica heterocedasticidad.
g7 <- ggplot(datos_d, aes(x = ajustados, y = residuos)) +
  geom_point(color = "#2980B9", size = 2.5, alpha = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "loess", se = FALSE, color = "orange", linewidth = 0.9) +
  labs(title    = "Gráfica 7. Residuos vs Valores Ajustados",
       subtitle = "Patrón sistemático = heterocedasticidad o mala especificación",
       x = "Valores Ajustados", y = "Residuos") +
  theme_minimal(base_size = 12)

ggsave(paste0(ruta_salida, "Grafica7_Residuos_vs_Ajustados.png"),
       g7, width = 8, height = 5, dpi = 150)

## `geom_smooth()` using formula = 'y ~ x'

print(g7)

## `geom_smooth()` using formula = 'y ~ x'

cat("Gráfica 7 exportada.\n\n")

## Gráfica 7 exportada.

# ------------------------------------------------------------------
# SUPUESTO 3 – NO AUTOCORRELACIÓN DE RESIDUOS
#
# Los errores de un período no deben estar correlacionados con los
# del período anterior. En series de tiempo esto es muy frecuente.
# Si hay autocorrelación, los errores estándar son ineficientes.
#
# Durbin-Watson:
#   DW ≈ 2 → sin autocorrelación ✓
#   DW < 2 → autocorrelación positiva (residuos del mismo signo se repiten)
#   DW > 2 → autocorrelación negativa
#
# Breusch-Godfrey: más potente, detecta autocorrelación de orden > 1
#   H0: no hay autocorrelación | p < 0.05 → hay autocorrelación
# ------------------------------------------------------------------

cat("--- Supuesto 3: Autocorrelación de residuos ---\n")

## --- Supuesto 3: Autocorrelación de residuos ---

dw <- dwtest(modelo1)
bg <- bgtest(modelo1, order = 2)

cat(paste("  Durbin-Watson:      DW =", round(dw$statistic, 4),
          "| p-valor =", round(dw$p.value, 4), "|",
          ifelse(dw$p.value > 0.05, "Sin autocorrelación ✓", "Hay autocorrelación"), "\n"))

##   Durbin-Watson:      DW = 1.5696 | p-valor = 0.1456 | Sin autocorrelación ✓

cat(paste("  Breusch-Godfrey(2): LM =", round(bg$statistic, 4),
          "| p-valor =", round(bg$p.value, 4), "|",
          ifelse(bg$p.value > 0.05, "Sin autocorrelación ✓", "Hay autocorrelación"), "\n\n"))

##   Breusch-Godfrey(2): LM = 2.2119 | p-valor = 0.3309 | Sin autocorrelación ✓

# Gráfica de residuos en el tiempo:
# No debe haber patrón cíclico (indicaría autocorrelación)
g8 <- ggplot(datos_d, aes(x = fecha, y = residuos)) +
  geom_line(color = "#8E44AD", linewidth = 0.9) +
  geom_point(size = 1.5, color = "#8E44AD") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title    = "Gráfica 8. Residuos en el Tiempo",
       subtitle = "Un patrón cíclico indicaría autocorrelación serial",
       x = NULL, y = "Residuos") +
  theme_minimal(base_size = 12)

ggsave(paste0(ruta_salida, "Grafica8_Residuos_Tiempo.png"),
       g8, width = 10, height = 5, dpi = 150)
print(g8)

cat("Gráfica 8 exportada.\n\n")

## Gráfica 8 exportada.

# ------------------------------------------------------------------
# SUPUESTO 4 – NO MULTICOLINEALIDAD (VIF)
#
# Las variables independientes no deben estar muy correlacionadas
# entre sí. Si lo están, el modelo no puede separar sus efectos
# individuales → coeficientes inestables con SE muy grandes.
#
# VIF (Factor de Inflación de la Varianza):
#   VIF < 5    → sin problema ✓
#   VIF 5–10   → multicolinealidad moderada (vigilar)
#   VIF > 10   → multicolinealidad grave → eliminar o combinar variables
# ------------------------------------------------------------------

cat("--- Supuesto 4: Multicolinealidad (VIF) ---\n")

## --- Supuesto 4: Multicolinealidad (VIF) ---

cat("VIF < 5: sin problema | 5-10: moderada | > 10: grave\n")

## VIF < 5: sin problema | 5-10: moderada | > 10: grave

vif_vals <- vif(modelo1)
print(round(vif_vals, 3))

##  d_homicidios   d_hurto_per   d_hurto_com   d_subempleo   d_ocupacion 
##         1.624         2.273         2.254         1.315         1.726 
##         d_PET d_defunciones d_nacimientos         d_IPC 
##         1.023         1.704         1.305         1.585

cat(paste("\n  VIF máximo:", round(max(vif_vals), 2),
          ifelse(max(vif_vals) < 5, "→ Sin problema ✓",
                 ifelse(max(vif_vals) < 10, "→ Moderada, vigilar",
                        "→ Grave, considerar eliminar variables")), "\n\n"))

## 
##   VIF máximo: 2.27 → Sin problema ✓

SECCIÓN 8: RESUMEN EJECUTIVO

# Imprime un resumen limpio con los resultados más importantes.
# ==============================================================================

s1 <- summary(modelo1)
s2 <- summary(modelo2)

cat("\n")

cat("================================================================\n")

## ================================================================

cat("               RESUMEN FINAL DE RESULTADOS\n")

##                RESUMEN FINAL DE RESULTADOS

cat("================================================================\n\n")

## ================================================================

cat("BASE DE DATOS:\n")

## BASE DE DATOS:

cat(paste("  Período:              2015 Q1 – 2024 Q3\n"))

##   Período:              2015 Q1 – 2024 Q3

cat(paste("  Obs. en niveles:      39\n"))

##   Obs. en niveles:      39

cat(paste("  Obs. en diferencias:  38\n\n"))

##   Obs. en diferencias:  38

cat("MODELO COMPLETO (modelo1) – 9 regresores:\n")

## MODELO COMPLETO (modelo1) – 9 regresores:

cat(paste("  R² ajustado:     ", round(s1$adj.r.squared, 4), "\n"))

##   R² ajustado:      0.8465

cat(paste("  F-estadístico:   ", round(s1$fstatistic[1], 3),
          "  p-valor:",
          round(pf(s1$fstatistic[1], s1$fstatistic[2],
                   s1$fstatistic[3], lower.tail = FALSE), 4), "\n\n"))

##   F-estadístico:    23.667   p-valor: 0

cat("MODELO PARSIMONIOSO (modelo2) – 5 regresores:\n")

## MODELO PARSIMONIOSO (modelo2) – 5 regresores:

cat(paste("  R² ajustado:     ", round(s2$adj.r.squared, 4), "\n"))

##   R² ajustado:      0.8184

cat(paste("  F-estadístico:   ", round(s2$fstatistic[1], 3),
          "  p-valor:",
          round(pf(s2$fstatistic[1], s2$fstatistic[2],
                   s2$fstatistic[3], lower.tail = FALSE), 4), "\n\n"))

##   F-estadístico:    34.344   p-valor: 0

cat("ARCHIVOS GENERADOS:\n")

## ARCHIVOS GENERADOS:

cat("  Tablas:\n")

##   Tablas:

cat("    Tabla1_Descriptivas.txt\n")

##     Tabla1_Descriptivas.txt

cat("    Tabla2_ADF_Niveles.csv\n")

##     Tabla2_ADF_Niveles.csv

cat("    Tabla3_ADF_Diferencias.csv\n")

##     Tabla3_ADF_Diferencias.csv

cat("    Tabla4_Modelo_Completo.txt\n")

##     Tabla4_Modelo_Completo.txt

cat("    Tabla5_Comparacion_Modelos.txt\n")

##     Tabla5_Comparacion_Modelos.txt

cat("  Gráficas:\n")

##   Gráficas:

cat("    Grafica1_Desempleo.png\n")

##     Grafica1_Desempleo.png

cat("    Grafica2_Seguridad.png\n")

##     Grafica2_Seguridad.png

cat("    Grafica3_MercadoLaboral.png\n")

##     Grafica3_MercadoLaboral.png

cat("    Grafica4_Diferenciacion.png\n")

##     Grafica4_Diferenciacion.png

cat("    Grafica5_Histograma_Residuos.png\n")

##     Grafica5_Histograma_Residuos.png

cat("    Grafica6_QQplot.png\n")

##     Grafica6_QQplot.png

cat("    Grafica7_Residuos_vs_Ajustados.png\n")

##     Grafica7_Residuos_vs_Ajustados.png

cat("    Grafica8_Residuos_Tiempo.png\n")

##     Grafica8_Residuos_Tiempo.png

cat("================================================================\n")

## ================================================================