#LIBRERIAS
#install.packages("readxl")
#install.packages("tidyverse")
#install.packages("lubridate")
#install.packages("car")      # Para pruebas de multicolinealidad
#install.packages("olsrr")  
# Cargar librerías
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
#library(olsrr)
#ARCHIVO
df <- read_excel("C:\\Users\\Diego Pérez\\Downloads\\Base_MTY_Final_Completo.xlsx")
# Verifica el contenido
glimpse(df)
## Rows: 48
## Columns: 18
## $ Fecha                   <dttm> 2020-01-01, 2020-02-01, 2020-03-01, 2020-04-0…
## $ Llegadas                <dbl> 5835, 5206, 4798, 2334, 2363, 2568, 3360, 3775…
## $ Salidas                 <dbl> 4167, 3810, 3156, 536, 401, 681, 1455, 1802, 2…
## $ Evento_Categoria        <dbl> 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1…
## $ PIB_Nuevo_Leon          <dbl> 145263.8, 145263.8, 145263.8, 145263.8, 145263…
## $ Poblacion_Nuevo_Leon    <dbl> 482036.8, 482036.8, 482036.8, 482036.8, 482036…
## $ Tipo_Cambio_MXN_USD     <dbl> 21.5, 21.5, 21.5, 21.5, 21.5, 21.5, 21.5, 21.5…
## $ IPC                     <dbl> 100.0000, 100.5319, 101.0638, 101.5957, 102.12…
## $ Movilidad_Urbana        <dbl> 3000.000, 3045.455, 3090.909, 3136.364, 3181.8…
## $ Ocupacion_Hotelera      <dbl> 30.00000, 31.36364, 32.72727, 34.09091, 35.454…
## $ INPC                    <dbl> 105.2, 105.4, 105.9, 106.3, 106.5, 106.8, 107.…
## $ Transporte_Urbano       <dbl> 2900, 2950, 3000, 2800, 2700, 2750, 3100, 3300…
## $ Ocupacion_Hotelera_REAL <dbl> 91249.02, 108145.57, 159893.79, 180163.90, 206…
## $ INPC_REAL               <dbl> 106.8890, 106.8890, 106.8380, 106.5000, 106.16…
## $ Evento_Pequeño          <dbl> 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1…
## $ Evento_Mediano          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ Evento_Magnoevento      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ Tipo_Evento             <chr> "Cultural", "Musical", "Deportivo", "Ninguno",…
df <- df %>%
  mutate(Fecha = as_date(Fecha))  # Asegura que sea clase Date
# Llegadas por año
df %>%
  mutate(Año = year(Fecha)) %>%
  group_by(Año) %>%
  summarise(Total_Llegadas = sum(Llegadas, na.rm = TRUE))
## # A tibble: 4 × 2
##     Año Total_Llegadas
##   <dbl>          <dbl>
## 1  2020          47061
## 2  2021          68542
## 3  2022          93340
## 4  2023         123764
# Gráfico de llegadas mensuales
df %>%
  ggplot(aes(x = Fecha, y = Llegadas)) +
  geom_line(color = "steelblue") +
  labs(title = "Llegadas aéreas a MTY por mes", x = "Fecha", y = "Llegadas")

#MODELO REGRESIÓN
# Eliminar columnas no numéricas o irrelevantes para regresión
df_model <- df %>%
  select(Fecha, Llegadas, Salidas, PIB_Nuevo_Leon, INPC_REAL, Ocupacion_Hotelera_REAL, Evento_Categoria, 
         Transporte_Urbano, Evento_Pequeño, Evento_Mediano, Evento_Magnoevento) %>%
  na.omit()
modelo_ols <- lm(Ocupacion_Hotelera_REAL ~ Llegadas + Salidas + PIB_Nuevo_Leon + INPC_REAL + Transporte_Urbano + Evento_Categoria, data = df)

summary(modelo_ols)
## 
## Call:
## lm(formula = Ocupacion_Hotelera_REAL ~ Llegadas + Salidas + PIB_Nuevo_Leon + 
##     INPC_REAL + Transporte_Urbano + Evento_Categoria, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -43624 -12315  -1038  12424  42141 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       564262.353 238663.949   2.364 0.022884 *  
## Llegadas              -3.050      4.842  -0.630 0.532242    
## Salidas              -22.296      7.011  -3.180 0.002800 ** 
## PIB_Nuevo_Leon         2.039      1.418   1.437 0.158187    
## INPC_REAL          -9593.026   2359.897  -4.065 0.000212 ***
## Transporte_Urbano    138.256     24.642   5.611 1.55e-06 ***
## Evento_Categoria      27.552   5149.730   0.005 0.995757    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19340 on 41 degrees of freedom
## Multiple R-squared:  0.6969, Adjusted R-squared:  0.6525 
## F-statistic: 15.71 on 6 and 41 DF,  p-value: 2.888e-09
vif(modelo_ols)
##          Llegadas           Salidas    PIB_Nuevo_Leon         INPC_REAL 
##         19.725455          5.930389         15.462199         51.676401 
## Transporte_Urbano  Evento_Categoria 
##         45.664068          1.820666
modelo_limpio <- lm(Ocupacion_Hotelera_REAL ~ Salidas + Transporte_Urbano, data = df_model)
summary(modelo_limpio)
## 
## Call:
## lm(formula = Ocupacion_Hotelera_REAL ~ Salidas + Transporte_Urbano, 
##     data = df_model)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -69107 -13943  -1824  17340  48163 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       120897.746  19904.189   6.074 2.42e-07 ***
## Salidas              -18.707      5.365  -3.487   0.0011 ** 
## Transporte_Urbano     40.486      6.796   5.957 3.62e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24840 on 45 degrees of freedom
## Multiple R-squared:  0.4508, Adjusted R-squared:  0.4263 
## F-statistic: 18.47 on 2 and 45 DF,  p-value: 1.395e-06

DATOS PANEL`

# Librerías necesarias
#install.packages("plm")
library(readxl)
library(plm)
## Warning: package 'plm' was built under R version 4.3.3
## 
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead
library(dplyr)

# Cargar el archivo
df <- read_excel("C:\\Users\\Diego Pérez\\Downloads\\Base_MTY_Final_Completo.xlsx")

# Convertir Fecha a tipo Date
df$Fecha <- as.Date(df$Fecha)

# Asegurar que Evento_Categoria esté como factor (ID del panel)
df$Categoria_Evento <- as.factor(df$Evento_Categoria)

# Crear subconjunto solo con variables que vamos a usar
df_model <- df %>%
  select(Fecha, Categoria_Evento, Ocupacion_Hotelera_REAL, Salidas, Llegadas, Transporte_Urbano) %>%
  filter(!is.na(Ocupacion_Hotelera_REAL))  # aseguramos no tener NAs

# Crear el panel
pdata <- pdata.frame(df_model, index = c("Categoria_Evento", "Fecha"))

MODELO DE EFECTOS FIJOS

modelo_fe <- plm(Ocupacion_Hotelera_REAL ~ Llegadas + Salidas + Transporte_Urbano,
                 data = pdata, model = "within")
summary(modelo_fe)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = Ocupacion_Hotelera_REAL ~ Llegadas + Salidas + 
##     Transporte_Urbano, data = pdata, model = "within")
## 
## Unbalanced Panel: n = 3, T = 5-37, N = 48
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -52175.7 -11693.8   1207.1  12109.7  47519.6 
## 
## Coefficients:
##                   Estimate Std. Error t-value  Pr(>|t|)    
## Llegadas          -15.7779     4.4916 -3.5128  0.001075 ** 
## Salidas            -2.7485     7.4572 -0.3686  0.714301    
## Transporte_Urbano  79.0286    12.4884  6.3281 1.339e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4.9654e+10
## Residual Sum of Squares: 2.0977e+10
## R-Squared:      0.57753
## Adj. R-Squared: 0.52724
## F-statistic: 19.1385 on 3 and 42 DF, p-value: 5.6329e-08

COMPARAR EFECTOS ALEATORIOS

#modelo_re <- plm(Ocupacion_Hotelera_REAL ~ Salidas + Transporte_Urbano, data = pdata, model = "random")
#summary(modelo_re)
#Error in swar_Between_check(estm[[2L]], method) : 
 # model not estimable: 3 coefficient(s) (incl. intercept) to be estimated but only 3 individual(s) in data for the between mode

Prueba Hausman

#phtest(modelo_fe, modelo_re)

PREDICCIONES 2024 Y 2025

# Simular proyección mensual de 2024 (ajusta según tu criterio o tasa de crecimiento)
#futuro <- data.frame(
#  Fecha = seq(as.Date("2024-01-01"), as.Date("2024-12-01"), by = "1 month"),
#  Categoria_Evento = factor(rep(c(1, 2, 3), length.out = 12)),  # puedes variar esto
#  Salidas = mean(df_model$Salidas) * 1.05,  # ejemplo: +5%
#  Transporte_Urbano = mean(df_model$Transporte_Urbano) * 1.05)

# Convertir a panel
#futuro_panel <- pdata.frame(futuro, index = c("Categoria_Evento", "Fecha"))

# Predecir (usando coeficientes del modelo)
#predicciones <- predict(modelo_fe, newdata = futuro_panel)

# Añadir a tabla
#futuro$Ocupacion_Predicha <- predicciones

ARIMA

#install.packages("forecast")
#install.packages("tseries")
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.1
# Usar solo la serie temporal de ocupación
library(readxl)
df <- read_excel("C:\\Users\\Diego Pérez\\Downloads\\Base_MTY_Final_Completo.xlsx")
df$Fecha <- as.Date(df$Fecha)

# Ordenar por fecha por si acaso
df <- df[order(df$Fecha), ]

# Crear la serie temporal mensual
ts_ocupacion <- ts(df$Llegadas, start = c(2020, 1), frequency = 12)

# Visualizar
plot(ts_ocupacion, main="Llegadas MTY (2020–2023)", ylab="Ocupación", xlab="Año")

# Descomposición clásica (modelo aditivo)
descomposicion <- decompose(ts_ocupacion, type = "additive")

# Mostrar los componentes
plot(descomposicion)

# Descomposición STL (más flexible)
descomposicion_stl <- stl(ts_ocupacion, s.window = "periodic")
plot(descomposicion_stl, main = "Descomposición STL de Llegadas MTY")

modelo_arima <- auto.arima(ts_ocupacion)
summary(modelo_arima)
## Series: ts_ocupacion 
## ARIMA(0,1,1)(0,1,0)[12] 
## 
## Coefficients:
##           ma1
##       -0.6346
## s.e.   0.1540
## 
## sigma^2 = 665820:  log likelihood = -284.07
## AIC=572.13   AICc=572.51   BIC=575.24
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 164.3697 686.7476 444.2707 2.560648 6.342476 0.2049613 -0.07602481
prediccion <- forecast(modelo_arima, h = 24)

# Visualizar
plot(prediccion, main="Pronóstico de Llegadas en 2026", xlab="Meses", ylab="Llegadas Estimada")

pred_df <- data.frame(
  Mes = seq(as.Date("2024-01-01"), by = "month", length.out = 24),
  Llegadas_Pronosticada = as.numeric(prediccion$mean)
)
print(pred_df)
##           Mes Llegadas_Pronosticada
## 1  2024-01-01              11977.59
## 2  2024-02-01              13768.59
## 3  2024-03-01              12276.59
## 4  2024-04-01              13298.59
## 5  2024-05-01              11486.59
## 6  2024-06-01              11623.59
## 7  2024-07-01              11969.59
## 8  2024-08-01              13935.59
## 9  2024-09-01              14264.59
## 10 2024-10-01              13217.59
## 11 2024-11-01              14195.59
## 12 2024-12-01              14108.59
## 13 2025-01-01              14674.18
## 14 2025-02-01              16465.18
## 15 2025-03-01              14973.18
## 16 2025-04-01              15995.18
## 17 2025-05-01              14183.18
## 18 2025-06-01              14320.18
## 19 2025-07-01              14666.18
## 20 2025-08-01              16632.18
## 21 2025-09-01              16961.18
## 22 2025-10-01              15914.18
## 23 2025-11-01              16892.18
## 24 2025-12-01              16805.18

ARIMAX

# Serie temporal de llegadas
ts_llegadas <- ts(df$Llegadas, start = c(2020, 1), frequency = 12)

# Variable exógena: transporte urbano
xreg <- df$Transporte_Urbano

# Sumar total de eventos por mes
df$Total_Eventos <- df$Evento_Pequeño + df$Evento_Mediano + df$Evento_Magnoevento

# Crear matriz de regresores exógenos
xreg <- as.matrix(df[, c("Transporte_Urbano", "Total_Eventos")])
# Ajustar modelo ARIMA con regresor exógeno
modelo_arimax <- auto.arima(ts_llegadas, xreg = xreg)
summary(modelo_arimax)
## Series: ts_llegadas 
## Regression with ARIMA(1,0,0)(0,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sma1  intercept  Transporte_Urbano  Total_Eventos
##       0.5414  0.3477  -5189.861             2.9484       140.2154
## s.e.  0.1418  0.1746   1380.449             0.3379       398.3954
## 
## sigma^2 = 614124:  log likelihood = -386.29
## AIC=784.57   AICc=786.62   BIC=795.8
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -45.85309 741.7229 565.1814 -2.958775 9.285318 0.2607426
##                    ACF1
## Training set -0.1259129
# Proyectar transporte urbano a futuro (12 meses)
#ultimo_promedio <- mean(tail(xreg, 12))
#xreg_futuro <- rep(ultimo_promedio * 1.05, 12)  # proyección plana +5%

# Hacer predicción
#pred_arimax <- forecast(modelo_arimax, xreg = xreg_futuro, h = 12)

#plot(pred_arimax, main = "Pronóstico de Llegadas con Transporte Urbano como regresor (ARIMAX)",
#     ylab = "Llegadas Estimadas", xlab = "Meses")

ARIMAX CON EVENTOS

ts_llegadas <- ts(df$Llegadas, start = c(2020, 1), frequency = 12)

modelo_arimax2 <- auto.arima(ts_llegadas, xreg = xreg)
summary(modelo_arimax2)
## Series: ts_llegadas 
## Regression with ARIMA(1,0,0)(0,0,1)[12] errors 
## 
## Coefficients:
##          ar1    sma1  intercept  Transporte_Urbano  Total_Eventos
##       0.5414  0.3477  -5189.861             2.9484       140.2154
## s.e.  0.1418  0.1746   1380.449             0.3379       398.3954
## 
## sigma^2 = 614124:  log likelihood = -386.29
## AIC=784.57   AICc=786.62   BIC=795.8
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -45.85309 741.7229 565.1814 -2.958775 9.285318 0.2607426
##                    ACF1
## Training set -0.1259129
# Promedios actuales
prom_transporte <- mean(tail(df$Transporte_Urbano, 12))
prom_eventos <- mean(tail(df$Total_Eventos, 12))

# Matriz de proyecciones exógenas para 2024 (12 meses)
xreg_futuro <- matrix(c(rep(prom_transporte * 1.05, 12),
                        rep(prom_eventos, 12)), ncol = 2)

pred_final <- forecast(modelo_arimax2, xreg = xreg_futuro, h = 12)
## Warning in forecast.forecast_ARIMA(modelo_arimax2, xreg = xreg_futuro, h = 12):
## xreg contains different column names from the xreg used in training. Please
## check that the regressors are in the same order.
plot(pred_final,
     main = "Pronóstico de Llegadas a MTY (ARIMAX con Transporte y Eventos)",
     xlab = "Meses", ylab = "Llegadas estimadas")