Limpieza de datos

Importamos la base de datos y la almacenamos como una dataframe. Posteriormente, aprovechamos la paquetería de dplyr para realizar la extracción de los datos que nos servirán en el análisis. Sacamos columnas como defunciones y viviendas dañadas y quitamos capturas con registros vacíos o no no válidos.

Analizamos todos los eventos dentro de una sola clasificación. Al final planteamos la posibilidad de segmentar los datos por estado con el fin de homogenizar.

Basehistorica_2000_a_2023 <- read_excel("Basehistorica_ 2000_a_2023.xlsx")
df <- Basehistorica_2000_a_2023
df <- dplyr::select(df, -c("Fecha de Inicio", "Fecha de Fin", "Descripcion general de los daños", 
                           "Fuente", "Población afectada", "Defunciones", "Viviendas dañadas", 
                           "Escuelas", "Hospitales", "Comercios", "Area de cultivo dañada / pastizales (h)"))

df_filtrado <- df[df$'Clasificación del fenómeno' == "Hidrometeorológico", ]

df_filtrado <- df_filtrado %>%
  filter(is.finite(`Total de daños (millones de pesos)`) & 
           !is.na(`Total de daños (millones de pesos)`) &
           is.numeric(`Total de daños (millones de pesos)`))

Supuestos del modelo

Para homogeneizar los valores de pérdidas económicas a lo largo del tiempo, se aplicó un ajuste inflacionario. Se asumió que todos los eventos ocurrieron al inicio del año de su registro y se actualizaron a valor presente de 2023 utilizando una tasa acumulada de inflación.

Además, se eliminaron eventos con pérdidas menores a 1 millón de pesos. Esto con la intención de evitar los valores negativos al aplicarle logaritmo natural a las pérdidas.

tasa_inflacion_acum <- c(1.0466, 1.12844412, 1.211497607, 1.249659782, 1.285025154, 1.347091869, 
                          1.438289988, 1.486616532, 1.518281464, 1.580227348, 1.642962373, 1.70161613, 
                          1.766617866, 1.844349052, 1.910192313, 2.034927871, 2.111441159, 2.196954526, 
                          2.270113112, 2.387931983, 2.482971676, 2.624501061, 2.739979108, 2.985481236)
 
 # Aplicar factores de crecimiento
 df_filtrado <- df_filtrado %>% 
   mutate(
     Factor_crecimiento = tasa_inflacion_acum[2024 - Año],
     `Daños proyectados a 2023` = `Total de daños (millones de pesos)` * Factor_crecimiento,
     `Log daños proyectados` = log(`Daños proyectados a 2023`)
   ) %>%
   filter(`Daños proyectados a 2023` > 1)

Análisis descriptivo de los datos

Antes de realizar los ajustes de modelos probabilísticos, se presentan algunas estadísticas descriptivas clave. Se calculan medidas de tendencia central, dispersión y forma para evaluar la distribución de los datos.

## Mínimo: 1.001476
## Máximo: 88131.07
## Rango: 88130.07
## Cuartiles:
##        25%        50%        75% 
##   4.097064  18.309377 185.400534
## Mediana: 18.30938
## Rango intercuartílico: 181.3035
## Valor en el percentil 99 (VaR.a): 9406.471
## Media: 635.9769
## Varianza: 16066511
## Desviación estándar: 4008.305
## Coeficiente de variación: 6.302596
## Skewness: 14.11005
## Kurtosis: 242.3429

Se presentan también visualizaciones de los datos, como histogramas, diagramas de caja y funciones de distribución empírica.

Modelos

Para encontrar una distribución que represente adecuadamente los datos de pérdidas, se ajustan cuatro modelos probabilísticos:

Distribución Normal

Distribución Log-Normal

Distribución Gamma

Distribución Weibull

Cada modelo es evaluado visualmente mediante gráficos de densidad, funciones de distribución acumulada y pruebas de bondad de ajuste.

Normal

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## mean 3.381508 0.06156649
## sd   2.377296 0.04353405
## Loglikelihood:  -3406.789   AIC:  6817.578   BIC:  6828.193 
## Correlation matrix:
##      mean sd
## mean    1  0
## sd      0  1

Log normal

## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters : 
##          estimate Std. Error
## meanlog 0.8261708 0.02864618
## sdlog   1.1061286 0.02025584
## Loglikelihood:  -3497.849   AIC:  6999.699   BIC:  7010.313 
## Correlation matrix:
##         meanlog sdlog
## meanlog       1     0
## sdlog         0     1

Gamma

## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##        estimate Std. Error
## shape 1.4186487 0.04709983
## rate  0.4196039 0.01665259
## Loglikelihood:  -3257.023   AIC:  6518.046   BIC:  6528.66 
## Correlation matrix:
##           shape      rate
## shape 1.0000000 0.8365621
## rate  0.8365621 1.0000000

Weibull

## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##       estimate Std. Error
## shape 1.323809 0.02815448
## scale 3.652218 0.07483484
## Loglikelihood:  -3230.608   AIC:  6465.217   BIC:  6475.831 
## Correlation matrix:
##           shape     scale
## shape 1.0000000 0.2974113
## scale 0.2974113 1.0000000

Resumen de todos los modelos

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.001475  1.410244  2.907413  3.381508  5.222518 11.386580

Validación de los modelos

Aquí vemos cómo se comparan los modelos con los datos registrados

summary(df_severidad)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.001475  1.410244  2.907413  3.381508  5.222518 11.386580
hist(df_severidad, breaks = 30, col = "lightblue", probability = TRUE, 
     main = "Histograma de severidad", xlab = "Severidad")
lines(density(df_severidad), col = "red", lwd = 2)

boxplot(df_severidad, main = "Boxplot de Severidad", col = "lightgreen")

plot.legend <- c("Normal", "Lognormal", "Gamma", "Weibull")
denscomp(list(ajuste_norm, ajuste_lognorm, ajuste_gamma, ajuste_weibull),
         legendtext = plot.legend)

qqcomp(list(ajuste_norm, ajuste_lognorm, ajuste_gamma, ajuste_weibull), 
       legendtext = c("Normal", "Lognormal", "Gamma", "Weibull"))

Aquí tenemos los resultados de la prueba de bondad entre los modelos

Anderson Darling Test de la normal:

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 3.38150809790743, sd = 2.37729575778577
##  Parameters assumed to be fixed
## 
## data:  df_severidad
## An = 22.629, p-value = 4.024e-07

Anderson Darling Test de la log normal:

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  with parameters meanlog = 0.826170791497832, sdlog = 1.10612856260064
##  Parameters assumed to be fixed
## 
## data:  df_severidad
## An = 39.85, p-value = 4.024e-07

Anderson Darling Test de la gamma:

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  with parameters shape = 1.41864874403211, rate = 0.419603931653432
##  Parameters assumed to be fixed
## 
## data:  df_severidad
## An = 11.655, p-value = 4.825e-07

Anderson Darling Test de la Weibull:

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.32380862106999, scale = 3.65221819203353
##  Parameters assumed to be fixed
## 
## data:  df_severidad
## An = 8.3146, p-value = 7.784e-05

Modelo elegido

La distribución de Weibull fue la que nos ofreció el valor más bajo tanto en el criterio de Akaike como en el Bayesiano, por lo que es el que eligiremos para modelar estos eventos.

Determinar el promedio por evento

Modelo General

Esta es la pérdida esperada si ocurre un evento hidrometereológico (en millones de pesos):

##       shape       scale 
##   0.3927605 102.0503477
## [1] 357.0978
## [1] 380.4611

Modelo considerando estados

resultados_estado <- df_filtrado %>%
  group_by(Estado) %>%
  summarise(
    Media = mean(`Log daños proyectados`),
    Mediana = median(`Log daños proyectados`),
    Varianza = var(`Log daños proyectados`)
  )

# Ver resultados
print(resultados_estado)
## # A tibble: 35 × 4
##    Estado                                       Media Mediana Varianza
##    <chr>                                        <dbl>   <dbl>    <dbl>
##  1 Aguascalientes                                2.60    2.41     5.49
##  2 Baja California                               4.04    3.50     6.18
##  3 Baja California Sur                           4.36    3.75     5.43
##  4 Campeche                                      4.46    4.65     6.19
##  5 Chiapas                                       3.05    2.39     6.32
##  6 Chiapas, Michoacan, Nuevo Leon, Quintana Roo  6.99    6.99    NA   
##  7 Chihuahua                                     3.66    3.62     3.91
##  8 Ciudad de México                              2.59    2.11     4.33
##  9 Coahuila                                      3.37    2.89     3.66
## 10 Colima                                        4.37    5.08     6.80
## # ℹ 25 more rows
# Visualización por estado
ggplot(df_filtrado, aes(x = `Log daños proyectados`, fill = Estado)) +
  geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  labs(title = "Distribución de severidad por Estado",
       x = "Daños proyectados (millones de pesos)",
       y = "Frecuencia") +
  theme_minimal()