1 Planteamiento del problema

El Instituto Costarricense de Turismo (ICT) requiere anticipar la llegada mensual de turistas internacionales por vía aérea durante los próximos doce meses, con el fin de apoyar la planificación de capacidad hotelera y aeroportuaria. El turismo es una de las principales fuentes de divisas del país, y las decisiones de capacidad —contratación de personal de migración, asignación de puertas y slots aeroportuarios, apertura de inventario hotelero— se toman con varios meses de antelación, por lo que dependen críticamente de la calidad del pronóstico.

Pregunta de decisión gerencial. ¿Cuál es el pronóstico de llegadas para los próximos cuatro trimestres y qué recomendación de capacidad debería trasladarse a las autoridades turísticas?

Reto técnico específico. La serie presenta (a) una estacionalidad anual muy marcada entre temporada alta (diciembre–abril) y temporada baja (mayo–noviembre) y (b) un quiebre estructural severo asociado al cierre de fronteras durante 2020–2021. La decisión metodológica central de este trabajo es cómo tratar ese período atípico antes de ajustar los modelos.

2 Datos

2.1 Fuente y carga

La serie corresponde a las llegadas mensuales de turistas internacionales por vía aérea a Costa Rica, publicadas por el Instituto Costarricense de Turismo con base en registros de la Dirección General de Migración y Extranjería. El período cubre de enero de 2019 a diciembre de 2025 (84 observaciones mensuales, siete años completos).

leer_datos <- function(ruta) {
  if (!file.exists(ruta)) return(NULL)
  bruto <- readxl::read_excel(ruta, range = readxl::cell_cols("A:B"))
  names(bruto)[1:2] <- c("fecha", "llegadas")
  bruto |>
    dplyr::filter(!is.na(llegadas)) |>
    dplyr::mutate(
      fecha = if (inherits(fecha, "Date") || inherits(fecha, "POSIXct")) {
        lubridate::floor_date(as.Date(fecha), "month")
      } else {
        as.Date(paste0(substr(as.character(fecha), 1, 7), "-01"))
      },
      llegadas = as.numeric(llegadas)
    ) |>
    dplyr::arrange(fecha)
}

datos <- leer_datos(params$ruta_datos)

# Respaldo embebido: garantiza que el documento sea reproducible aun si cambia la ruta.
if (is.null(datos)) {
  llegadas_v <- c(
    248754, 250939, 276036, 198989, 164992, 196977, 213125, 174109, 112102, 135320, 191346, 255611,
    261658, 276970, 127201,    505,    517,    625,    418,   1636,   3371,   9888,  36044,  71000,
     60619,  54745,  88391,  85846,  90321, 117527, 146915, 114227,  77599,  91292, 141099, 201902,
    169585, 184358, 225411, 194181, 160944, 186097, 204878, 159753, 102705, 116925, 181721, 231402,
    240487, 250038, 266578, 210717, 168248, 202235, 223296, 174913, 111696, 129201, 212646, 281095,
    270712, 291152, 322400, 221573, 199154, 227452, 239368, 182800, 105100, 122188, 201152, 278437,
    266675, 270810, 312844, 231678, 189881, 217120, 242408, 182342, 104246, 129369, 225679, 316226
  )
  datos <- data.frame(
    fecha = seq(as.Date("2019-01-01"), by = "month", length.out = length(llegadas_v)),
    llegadas = llegadas_v
  )
}

y <- ts(datos$llegadas, start = c(2019, 1), frequency = 12)
c(n = length(y), inicio = start(y)[1], fin = end(y)[1])
##      n inicio    fin 
##     84   2019   2025

2.2 Exploración inicial

autoplot(y) +
  annotate("rect", xmin = 2020.17, xmax = 2022, ymin = -Inf, ymax = Inf,
           alpha = .12, fill = "red") +
  annotate("text", x = 2021.1, y = max(y) * .92,
           label = "Quiebre estructural\n(cierre de fronteras)", size = 3) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", title = NULL)
Figura 1. Llegadas mensuales de turistas internacionales por vía aérea, Costa Rica, 2019–2025.

Figura 1. Llegadas mensuales de turistas internacionales por vía aérea, Costa Rica, 2019–2025.

Nota. Elaboración propia con datos del ICT.

y_post <- window(y, start = c(2022, 1))
ggsubseriesplot(y_post) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Mes", y = "Llegadas (personas)", title = NULL)
Figura 2. Gráfico de subseries estacionales del período posterior al quiebre (2022–2025).

Figura 2. Gráfico de subseries estacionales del período posterior al quiebre (2022–2025).

Nota. La línea horizontal de cada panel corresponde a la media del mes.

datos |>
  dplyr::mutate(anio = lubridate::year(fecha)) |>
  dplyr::group_by(anio) |>
  dplyr::summarise(total = sum(llegadas), .groups = "drop") |>
  dplyr::mutate(variacion = scales::percent(total / dplyr::lag(total) - 1, accuracy = .1)) |>
  knitr::kable(caption = "Tabla 1. Llegadas anuales por vía aérea y variación interanual.",
               format.args = list(big.mark = ","), col.names = c("Año", "Total", "Variación"))
Tabla 1. Llegadas anuales por vía aérea y variación interanual.
Año Total Variación
2,019 2,418,300 NA
2,020 789,833 -67.3%
2,021 1,270,483 60.9%
2,022 2,117,960 66.7%
2,023 2,471,150 16.7%
2,024 2,661,488 7.7%
2,025 2,689,278 1.0% 
indices <- data.frame(mes = cycle(y_post), valor = as.numeric(y_post)) |>
  dplyr::group_by(mes) |>
  dplyr::summarise(indice = mean(valor) / mean(as.numeric(y_post)), .groups = "drop") |>
  dplyr::mutate(mes = month.name[mes])
knitr::kable(indices, digits = 3,
             caption = "Tabla 2. Índice estacional promedio 2022–2025 (media global = 1).",
             col.names = c("Mes", "Índice"))
Tabla 2. Índice estacional promedio 2022–2025 (media global = 1).
Mes Índice
January 1.144
February 1.203
March 1.361
April 1.036
May 0.867
June 1.006
July 1.099
August 0.845
September 0.512
October 0.601
November 0.991
December 1.337

La exploración confirma tres hechos relevantes para el modelado: (a) el mes de setiembre opera cerca de la mitad del promedio anual, mientras marzo y diciembre lo superan en más de un tercio; (b) la amplitud de la oscilación estacional crece con el nivel de la serie, lo cual anticipa una estructura multiplicativa; y (c) el período comprendido entre marzo de 2020 y diciembre de 2021 no responde al proceso generador habitual de la serie.

3 Tratamiento del quiebre estructural 2020–2021

Ajustar cualquier modelo sobre la serie sin corregir el período de cierre de fronteras produce estimaciones sesgadas: el episodio actúa como un valor atípico extremo que infla la varianza residual y distorsiona los coeficientes estacionales. Se evaluaron tres estrategias.

Estrategia Descripción Ventaja Desventaja
A. Submuestra Utilizar solo 2022–2025 No requiere supuestos sobre el período atípico Solo 48 observaciones, 4 ciclos estacionales
B. Imputación Marcar 2020-03 a 2021-12 como faltante e imputar con interpolación estacional (na.interp) Conserva 84 observaciones y la estacionalidad de 2019 Introduce valores contrafactuales no observados
C. Variables de intervención Serie completa con variables dicotómicas COVID en xreg Modela explícitamente el choque Requiere supuestos sobre la duración del efecto

Se adopta la estrategia B como especificación principal y la estrategia A como verificación de robustez, comparando ambas con las mismas métricas fuera de muestra.

fechas <- datos$fecha
covid <- fechas >= as.Date(params$inicio_covid) & fechas <= as.Date(params$fin_covid)

y_na  <- y; y_na[covid] <- NA
y_imp <- forecast::na.interp(y_na)   # Estrategia B
y_sub <- window(y, start = c(2022, 1))  # Estrategia A

cat("Meses tratados como atípicos:", sum(covid), "\n")
## Meses tratados como atípicos: 22
autoplot(y_imp, series = "Imputada (Estrategia B)") +
  autolayer(y, series = "Observada") +
  scale_colour_manual(values = c("Observada" = "grey30", "Imputada (Estrategia B)" = "#2E7F79")) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", colour = NULL, title = NULL) +
  theme(legend.position = "bottom")
Figura 3. Serie observada frente a la serie con el período 2020–2021 imputado.

Figura 3. Serie observada frente a la serie con el período 2020–2021 imputado.

4 Descomposición

Dado que la amplitud estacional se expande junto con el nivel de la serie, se aplica una descomposición multiplicativa. De manera equivalente, se emplea la descomposición STL sobre el logaritmo de la serie, lo que permite una tendencia flexible y componentes estacionales que evolucionan lentamente.

desc_mult <- decompose(y_imp, type = "multiplicative")
autoplot(desc_mult) + labs(x = "Año", title = NULL)
Figura 4. Descomposición multiplicativa clásica de la serie imputada.

Figura 4. Descomposición multiplicativa clásica de la serie imputada. 

stl_log <- stl(log(y_imp), s.window = "periodic", robust = TRUE)
autoplot(stl_log) + labs(x = "Año", title = NULL)
Figura 5. Descomposición STL sobre el logaritmo de la serie imputada.

Figura 5. Descomposición STL sobre el logaritmo de la serie imputada. 

res_mult <- na.omit(desc_mult$random)
res_adit <- na.omit(decompose(y_imp, type = "additive")$random)
data.frame(
  Descomposicion = c("Multiplicativa", "Aditiva"),
  Var_residuo_relativa = c(sd(res_mult - 1), sd(res_adit) / mean(y_imp))
) |>
  knitr::kable(digits = 4,
    caption = "Tabla 3. Dispersión relativa del residuo según tipo de descomposición.")
Tabla 3. Dispersión relativa del residuo según tipo de descomposición.
Descomposicion Var_residuo_relativa
Multiplicativa 0.096
Aditiva 0.045

La descomposición multiplicativa deja un residuo con dispersión relativa menor y sin patrón sistemático por nivel, lo que justifica su elección frente a la aditiva. En consecuencia, todos los modelos posteriores se estiman con transformación logarítmica (lambda = 0), equivalente a un tratamiento multiplicativo de la estacionalidad.

5 Suavizamiento

mm12 <- forecast::ma(y_imp, order = 12, centre = TRUE)
autoplot(y_imp, series = "Serie imputada") +
  autolayer(mm12, series = "Media móvil (12)", size = 1) +
  scale_colour_manual(values = c("Serie imputada" = "grey60", "Media móvil (12)" = "#C2871F")) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", colour = NULL, title = NULL) +
  theme(legend.position = "bottom")
Figura 6. Serie original y media móvil centrada de orden 12.

Figura 6. Serie original y media móvil centrada de orden 12.

La media móvil de orden 12 elimina por completo la oscilación estacional y revela la trayectoria de recuperación: caída, rebote acelerado durante 2022–2023 y estabilización en un nivel superior al previo a la pandemia a partir de 2024.

hw_mult <- hw(y_imp, seasonal = "multiplicative", h = params$h)
ses_simple <- ses(y_imp, h = params$h)

data.frame(
  Metodo = c("Suavizamiento exponencial simple (SES)", "Holt-Winters multiplicativo"),
  AICc = c(ses_simple$model$aicc, hw_mult$model$aicc),
  RMSE_ajuste = c(accuracy(ses_simple)[, "RMSE"], accuracy(hw_mult)[, "RMSE"])
) |>
  knitr::kable(digits = 1,
    caption = "Tabla 4. Comparación entre suavizamiento exponencial simple y Holt-Winters.",
    format.args = list(big.mark = ","))
Tabla 4. Comparación entre suavizamiento exponencial simple y Holt-Winters.
Metodo AICc RMSE_ajuste
Suavizamiento exponencial simple (SES) 2,188.9 47,860.0
Holt-Winters multiplicativo 2,071.0 11,591.4 
autoplot(y_imp, series = "Observada") +
  autolayer(fitted(hw_mult), series = "Holt-Winters") +
  scale_colour_manual(values = c("Observada" = "grey40", "Holt-Winters" = "#7C3AED")) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", colour = NULL, title = NULL) +
  theme(legend.position = "bottom")
Figura 7. Ajuste de Holt-Winters multiplicativo sobre la serie imputada.

Figura 7. Ajuste de Holt-Winters multiplicativo sobre la serie imputada.

El suavizamiento exponencial simple es incapaz de reproducir la oscilación anual —su pronóstico es una línea plana—, mientras que Holt-Winters multiplicativo la captura con fidelidad. Esto confirma que cualquier modelo candidato debe incorporar componente estacional, y anticipa el resultado de la comparación ARIMA frente a SARIMA.

6 Estacionariedad

ly <- log(y_imp)
adf <- tseries::adf.test(ly)
kpss <- tseries::kpss.test(ly)

data.frame(
  Prueba = c("Dickey-Fuller aumentado (H0: raíz unitaria)",
             "KPSS (H0: estacionariedad)"),
  Estadistico = c(adf$statistic, kpss$statistic),
  Valor_p = c(adf$p.value, kpss$p.value)
) |>
  knitr::kable(digits = 3, caption = "Tabla 5. Pruebas de estacionariedad sobre el logaritmo de la serie.")
Tabla 5. Pruebas de estacionariedad sobre el logaritmo de la serie.
Prueba Estadistico Valor_p
Dickey-Fuller Dickey-Fuller aumentado (H0: raíz unitaria) -4.694 0.01
KPSS Level KPSS (H0: estacionariedad) 0.186 0.10 
c(diferencias_regulares = forecast::ndiffs(ly),
  diferencias_estacionales = forecast::nsdiffs(ly))
##    diferencias_regulares diferencias_estacionales 
##                        0                        1
ggtsdisplay(diff(diff(ly), lag = 12), main = NULL)
Figura 8. Serie en logaritmos tras diferenciación regular y estacional, con sus funciones de autocorrelación.

Figura 8. Serie en logaritmos tras diferenciación regular y estacional, con sus funciones de autocorrelación.

Las pruebas y los criterios automáticos coinciden en la necesidad de una diferenciación estacional; los correlogramas de la serie doblemente diferenciada muestran un rezago significativo en el retardo 1 y otro en el retardo 12, patrón compatible con una estructura de medias móviles regular y estacional.

7 Modelos ARIMA y SARIMA

7.1 Diseño de la validación

Se reserva el año 2025 completo (12 observaciones) como conjunto de prueba. Los modelos se estiman únicamente con información hasta diciembre de 2024 y se comparan por su error de pronóstico fuera de muestra, criterio más exigente que el ajuste dentro de muestra.

corte <- c(2024, 12)
train_B <- window(y_imp, end = corte)
train_A <- window(y_sub, end = corte)
test    <- window(y,     start = c(2025, 1))
c(train_B = length(train_B), train_A = length(train_A), test = length(test))
## train_B train_A    test 
##      72      36      12

7.2 Estimación

# --- Estrategia B (serie completa imputada) ---
arima_B  <- auto.arima(train_B, lambda = 0, seasonal = FALSE,
                       stepwise = FALSE, approximation = FALSE)
sarima_B <- auto.arima(train_B, lambda = 0, seasonal = TRUE,
                       stepwise = FALSE, approximation = FALSE)
manual_B <- Arima(train_B, order = c(0, 1, 1), seasonal = c(0, 1, 1), lambda = 0)

# --- Estrategia A (submuestra 2022-2024) ---
arima_A  <- auto.arima(train_A, lambda = 0, seasonal = FALSE)
sarima_A <- auto.arima(train_A, lambda = 0, seasonal = TRUE)

# --- Referencias de comparación ---
hw_train  <- hw(train_B, seasonal = "multiplicative", h = length(test))
snaive_tr <- snaive(train_B, h = length(test))
data.frame(
  Estrategia = c("B", "B", "B", "A", "A"),
  Modelo = c("ARIMA automático (sin estacionalidad)", "SARIMA automático",
             "SARIMA manual (0,1,1)(0,1,1)[12]",
             "ARIMA automático (sin estacionalidad)", "SARIMA automático"),
  Especificacion = c(as.character(arima_B), as.character(sarima_B),
                     as.character(manual_B), as.character(arima_A),
                     as.character(sarima_A))
) |>
  knitr::kable(caption = "Tabla 6. Especificaciones seleccionadas por estrategia.")
Tabla 6. Especificaciones seleccionadas por estrategia.
Estrategia Modelo Especificacion
B ARIMA automático (sin estacionalidad) ARIMA(3,0,2) with non-zero mean
B SARIMA automático ARIMA(2,1,0)(2,1,0)[12]
B SARIMA manual (0,1,1)(0,1,1)[12] ARIMA(0,1,1)(0,1,1)[12]
A ARIMA automático (sin estacionalidad) ARIMA(0,0,2) with non-zero mean
A SARIMA automático ARIMA(0,1,0)(0,1,0)[12]

7.3 Comparación de desempeño

evaluar <- function(objeto, etiqueta, prueba = test, h = length(test)) {
  pr <- if (inherits(objeto, "forecast")) objeto else forecast(objeto, h = h)
  a  <- forecast::accuracy(pr, prueba)
  data.frame(
    Modelo = etiqueta,
    AICc = ifelse(is.null(pr$model$aicc), NA, pr$model$aicc),
    RMSE = a["Test set", "RMSE"],
    MAE  = a["Test set", "MAE"],
    MAPE = a["Test set", "MAPE"]
  )
}

comparacion <- dplyr::bind_rows(
  evaluar(snaive_tr, "Naive estacional (referencia)"),
  evaluar(hw_train,  "Holt-Winters multiplicativo (B)"),
  evaluar(arima_B,   "ARIMA sin estacionalidad (B)"),
  evaluar(sarima_B,  "SARIMA automático (B)"),
  evaluar(manual_B,  "SARIMA manual (0,1,1)(0,1,1)[12] (B)"),
  evaluar(arima_A,   "ARIMA sin estacionalidad (A)"),
  evaluar(sarima_A,  "SARIMA automático (A)")
) |>
  dplyr::arrange(MAPE)

knitr::kable(comparacion, digits = c(0, 1, 0, 0, 2),
  caption = "Tabla 7. Desempeño fuera de muestra sobre el año 2025 (12 meses).",
  format.args = list(big.mark = ","))
Tabla 7. Desempeño fuera de muestra sobre el año 2025 (12 meses).
Modelo AICc RMSE MAE MAPE
Holt-Winters multiplicativo (B) 1,787.7 14,638 7,799 3.28
SARIMA manual (0,1,1)(0,1,1)[12] (B) -58.1 13,090 10,501 4.61
SARIMA automático (A) -58.4 15,958 11,270 4.67
Naive estacional (referencia) NA 15,566 11,458 4.73
SARIMA automático (B) -59.5 48,262 46,969 22.98
ARIMA sin estacionalidad (A) -16.7 63,775 51,177 25.60
ARIMA sin estacionalidad (B) -5.0 68,676 57,334 27.23 
mejor <- comparacion$Modelo[1]
cat("Modelo con menor MAPE fuera de muestra:", mejor, "\n")
## Modelo con menor MAPE fuera de muestra: Holt-Winters multiplicativo (B)
autoplot(window(y_imp, start = c(2023, 1))) +
  autolayer(forecast(sarima_B, h = 12), series = "SARIMA (B)", PI = FALSE, size = 1) +
  autolayer(forecast(arima_B,  h = 12), series = "ARIMA (B)",  PI = FALSE, size = 1) +
  autolayer(test, series = "Observado 2025", size = 1) +
  scale_colour_manual(values = c("SARIMA (B)" = "#2E7F79", "ARIMA (B)" = "#9B2C2C",
                                 "Observado 2025" = "black")) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", colour = NULL, title = NULL) +
  theme(legend.position = "bottom")
Figura 9. Pronóstico de 2025 según cada modelo frente a los valores observados.

Figura 9. Pronóstico de 2025 según cada modelo frente a los valores observados.

El contraste es contundente: el modelo sin componente estacional proyecta una trayectoria prácticamente plana y comete errores de magnitud comparable a la propia amplitud estacional, mientras que las especificaciones SARIMA reproducen tanto el pico de temporada alta como el valle de setiembre. El componente estacional no es un refinamiento marginal en esta serie, sino la fuente principal de su capacidad predictiva.

7.4 Diagnóstico de residuos

modelo_final_train <- sarima_B
checkresiduals(modelo_final_train)
Figura 10. Diagnóstico de residuos del modelo SARIMA seleccionado.

Figura 10. Diagnóstico de residuos del modelo SARIMA seleccionado. 

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)(2,1,0)[12]
## Q* = 4.0515, df = 10, p-value = 0.945
## 
## Model df: 4.   Total lags used: 14
lb <- Box.test(residuals(modelo_final_train), lag = 24, type = "Ljung-Box",
               fitdf = length(coef(modelo_final_train)))
sw <- shapiro.test(residuals(modelo_final_train))
data.frame(
  Prueba = c("Ljung-Box (24 rezagos)", "Shapiro-Wilk (normalidad)"),
  Estadistico = c(lb$statistic, sw$statistic),
  Valor_p = c(lb$p.value, sw$p.value)
) |>
  knitr::kable(digits = 4, caption = "Tabla 8. Pruebas sobre los residuos del modelo seleccionado.")
Tabla 8. Pruebas sobre los residuos del modelo seleccionado.
Prueba Estadistico Valor_p
X-squared Ljung-Box (24 rezagos) 11.6235 0.9284
W Shapiro-Wilk (normalidad) 0.8427 0.0000

Un valor p superior a 0.05 en la prueba de Ljung-Box indica que no queda autocorrelación sistemática en los residuos, es decir, que el modelo ha extraído la estructura temporal disponible. Debe reportarse el resultado obtenido y, si la normalidad se rechaza, advertir que los intervalos de predicción pueden estar levemente subestimados en las colas.

8 Pronóstico para 2026

El modelo seleccionado se reestima sobre la serie completa (2019–2025) y se proyectan doce meses con intervalos de predicción al 80 % y 95 %.

modelo_final <- auto.arima(y_imp, lambda = 0, seasonal = TRUE,
                           stepwise = FALSE, approximation = FALSE)
pron <- forecast(modelo_final, h = params$h, level = c(80, 95))
summary(modelo_final)
## Series: y_imp 
## ARIMA(2,1,0)(2,1,0)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1      ar2     sar1     sar2
##       0.1858  -0.3824  -0.6396  -0.5305
## s.e.  0.1092   0.1082   0.1040   0.1031
## 
## sigma^2 = 0.01364:  log likelihood = 48.56
## AIC=-87.12   AICc=-86.2   BIC=-75.81
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 3.516746 12893.67 9853.087 -0.3272685 6.508532 0.3298198 0.1376291
autoplot(pron) +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Año", y = "Llegadas (personas)", title = NULL)
Figura 11. Pronóstico de llegadas mensuales para 2026 con intervalos de predicción al 80 % y 95 %.

Figura 11. Pronóstico de llegadas mensuales para 2026 con intervalos de predicción al 80 % y 95 %. 

tabla_pron <- data.frame(
  Mes = format(seq(as.Date("2026-01-01"), by = "month", length.out = params$h), "%B %Y"),
  Pronostico = as.numeric(pron$mean),
  IC80_inf = as.numeric(pron$lower[, 1]), IC80_sup = as.numeric(pron$upper[, 1]),
  IC95_inf = as.numeric(pron$lower[, 2]), IC95_sup = as.numeric(pron$upper[, 2])
)
knitr::kable(tabla_pron, digits = 0,
  caption = "Tabla 9. Pronóstico mensual de llegadas para 2026.",
  format.args = list(big.mark = ","),
  col.names = c("Mes", "Pronóstico", "IC 80 % inf.", "IC 80 % sup.",
                "IC 95 % inf.", "IC 95 % sup."))
Tabla 9. Pronóstico mensual de llegadas para 2026.
Mes Pronóstico IC 80 % inf. IC 80 % sup. IC 95 % inf. IC 95 % sup.
January 2026 309,454 266,433 359,421 246,136 389,060
February 2026 311,459 246,925 392,860 218,366 444,241
March 2026 344,400 264,517 448,406 230,030 515,634
April 2026 264,839 199,350 351,842 171,518 408,934
May 2026 216,403 158,902 294,713 134,934 347,061
June 2026 253,006 180,824 354,002 151,368 422,889
July 2026 278,713 194,569 399,246 160,860 482,911
August 2026 214,859 146,956 314,137 120,188 384,102
September 2026 130,428 87,476 194,472 70,803 240,267
October 2026 154,765 101,814 235,253 81,571 293,635
November 2026 260,039 167,940 402,648 133,240 507,510
December 2026 352,842 223,917 556,000 176,011 707,329 
tabla_trim <- tabla_pron |>
  dplyr::mutate(Trimestre = rep(paste0("T", 1:4, " 2026"), each = 3)) |>
  dplyr::group_by(Trimestre) |>
  dplyr::summarise(Pronostico = sum(Pronostico),
                   IC80_inf = sum(IC80_inf), IC80_sup = sum(IC80_sup), .groups = "drop")

knitr::kable(tabla_trim, digits = 0,
  caption = "Tabla 10. Pronóstico trimestral de llegadas para 2026 (respuesta directa a la pregunta gerencial).",
  format.args = list(big.mark = ","),
  col.names = c("Trimestre", "Pronóstico", "IC 80 % inf.", "IC 80 % sup."))
Tabla 10. Pronóstico trimestral de llegadas para 2026 (respuesta directa a la pregunta gerencial).
Trimestre Pronóstico IC 80 % inf. IC 80 % sup.
T1 2026 965,313 777,875 1,200,687
T2 2026 734,248 539,075 1,000,556
T3 2026 624,000 429,001 907,855
T4 2026 767,646 493,670 1,193,900 
total_2026 <- sum(tabla_pron$Pronostico)
total_2025 <- sum(window(y, start = c(2025, 1)))
cat(sprintf("Total proyectado 2026: %s llegadas (%+.1f %% respecto de 2025: %s)\n",
            format(round(total_2026), big.mark = ","),
            100 * (total_2026 / total_2025 - 1),
            format(round(total_2025), big.mark = ",")))
## Total proyectado 2026: 3,091,207 llegadas (+14.9 % respecto de 2025: 2,689,278)

9 Recomendación gerencial

Redactar esta sección en lenguaje no técnico, respondiendo de forma directa la pregunta de decisión. Se sugiere estructurarla en cuatro puntos, sustituyendo los valores con los que arroje el modelo:

  1. Volumen esperado. El país recibiría alrededor de 3,091,207 turistas por vía aérea en 2026, con la mayor concentración en el primer trimestre.
  2. Meses de presión de capacidad. Enero, febrero, marzo y diciembre superan los 300 000 arribos mensuales proyectados; son los meses en que deben reforzarse puestos de migración, personal de rampa y disponibilidad de habitaciones.
  3. Meses de holgura. Setiembre y octubre operan cerca de la mitad del promedio anual, por lo que constituyen la ventana natural para mantenimiento de infraestructura, capacitación de personal y campañas de estímulo de demanda.
  4. Margen de planificación. La capacidad instalada debería dimensionarse contra el límite superior del intervalo al 80 % y no contra el valor central, dado que el costo de subestimar la demanda —congestión aeroportuaria, deterioro de la experiencia del visitante— es asimétricamente mayor que el de una holgura moderada.

10 Limitaciones del análisis

  • La imputación del período 2020–2021 genera valores contrafactuales que no fueron observados; aunque la comparación con la estrategia de submuestra respalda la decisión, la tasa de crecimiento proyectada es sensible a este supuesto.
  • El modelo es univariado: no incorpora variables explicativas como el tipo de cambio, la conectividad aérea, los precios de los competidores regionales ni las campañas de promoción del ICT.
  • La serie mide llegadas por vía aérea y excluye el ingreso terrestre y marítimo, por lo que el pronóstico no representa el total de visitantes al país.
  • Los intervalos de predicción suponen que la estructura estimada se mantiene estable; un choque externo de la magnitud del de 2020 quedaría por completo fuera de ellos.
  • El horizonte de doce meses implica una acumulación de incertidumbre creciente: la amplitud del intervalo del mes 12 es sustancialmente mayor que la del mes 1.

11 Referencias

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and practice (3.ª ed.). OTexts. https://otexts.com/fpp3/

Instituto Costarricense de Turismo. (2026). Anuario de estadísticas de turismo: Llegadas de turistas internacionales por vía aérea. https://www.ict.go.cr/

R Core Team. (2026). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0    knitr_1.51      urca_1.3-4      tseries_0.10-62
##  [5] forecast_9.0.2  lubridate_1.9.5 ggplot2_4.0.3   tidyr_1.3.2    
##  [9] dplyr_1.2.1     readxl_1.5.0   
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     lattice_0.22-9     digest_0.6.39     
##  [5] magrittr_2.0.5     evaluate_1.0.5     grid_4.6.1         timechange_0.4.0  
##  [9] RColorBrewer_1.1-3 fastmap_1.2.0      cellranger_1.1.0   jsonlite_2.0.0    
## [13] purrr_1.2.2        jquerylib_0.1.4    cli_3.6.6          rlang_1.3.0       
## [17] withr_3.0.3        cachem_1.1.0       yaml_2.3.12        tools_4.6.1       
## [21] parallel_4.6.1     colorspace_2.1-3   curl_7.1.0         vctrs_0.7.3       
## [25] R6_2.6.1           zoo_1.9-0          lifecycle_1.0.5    pkgconfig_2.0.3   
## [29] pillar_1.11.1      bslib_0.12.0       gtable_0.3.6       glue_1.8.1        
## [33] quantmod_0.4.29    Rcpp_1.1.2         xfun_0.60          tibble_3.3.1      
## [37] tidyselect_1.2.1   farver_2.1.2       htmltools_0.5.9    nlme_3.1-169      
## [41] labeling_0.4.3     rmarkdown_2.31     xts_0.14.2         timeDate_4052.112 
## [45] fracdiff_1.5-4     compiler_4.6.1     S7_0.2.2           quadprog_1.5-8    
## [49] TTR_0.24.4