Introducción

En este trabajo realizamos un análisis de series de tiempo aplicado al sector financiero, con el propósito de entender cómo se comportan en el tiempo cuatro variables económicas mensuales, identificar sus patrones de tendencia y estacionalidad, y construir un modelo capaz de pronosticar el comportamiento futuro de una de ellas. Para lograrlo, seguimos la metodología clásica de Box-Jenkins para modelos ARIMA/SARIMA.

#Cargar librerías necesarias
library(readxl)  # Para leer archivos Excel
library(tseries)  # Para pruebas de estacionariedad
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)  # Para modelado ARIMA y pronósticos
library(ggplot2)  # Para visualización de datos
library(plotly)  # Para gráficos interactivos
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(timetk)   #timetk simplifica y acelera el análisis exploratorio, visualización, y preparación de datos temporales para modelado. Es ideal para quienes trabajan con series temporales en un flujo de trabajo "tidy" y buscan integrar análisis visuales, detección de patrones y forecasting en un solo paquete.
## Warning: package 'timetk' was built under R version 4.4.3
library(readxl)
data_col <- read_excel("C:/Users/yijim/Downloads/Base_Caso.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))
View(data_col)

Variables

Declaramos cada una de estas cuatro series como objetos de serie de tiempo (ts) con periodicidad mensual, para que las herramientas de descomposición y modelado que usamos más adelante reconocieran correctamente la estructura temporal de los datos.

# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$TRM, start = c(2012, 1), frequency = 12)
# Convertir/declarar el ISE en serie de tiempo  mensual
variable2_ts <- ts(data_col$ISE, start = c(2012, 1), frequency = 12)
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(data_col$ICC, start = c(2012, 1), frequency = 12)
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable4_ts <- ts(data_col$POLLO, start = c(2012, 1), frequency = 12)

Extraccion señales

Para cada una de las cuatro variables, aplicamos una descomposición STL (Seasonal-Trend decomposition using Loess). Elegimos esta técnica porque separa la serie observada en tres componentes que nos interesaba analizar por separado: • Componente estacional: los patrones que se repiten todos los años en las mismas épocas. • Componente de tendencia: hacia dónde se mueve la variable en el mediano y largo plazo, una vez removemos el ruido de corto plazo. • Componente residual: lo que queda después de retirar tendencia y estacionalidad; idealmente, este componente debería comportarse como ruido sin patrones predecibles.

Serie original vs. tendencia

• La tendencia es claramente creciente y prácticamente monótona entre 2012 y 2020, se estanca brevemente durante la pandemia, retoma el crecimiento hasta 2023 y muestra una leve inflexión a la baja al cierre de 2024. • El componente estacional es sistemático y de magnitud considerable: oscila aproximadamente entre −9.000 y +7.000 unidades alrededor de la tendencia en cada ciclo anual, lo que confirma que existe un patrón de estacionalidad de calendario relevante en las exportaciones de pollo — algo que tuvimos en cuenta al momento de decidir si el modelo de pronóstico debía incluir un componente estacional. • El residuo permanece cercano a cero y sin patrones aparentes durante casi toda la muestra, con la notable excepción de 2020, donde alcanza un valor de casi −45.000 unidades: la caída de la pandemia fue un evento verdaderamente atípico, no explicado ni por la tendencia ni por la estacionalidad de la serie.

Variable 1 (original)

library(ggplot2)
library(plotly)

# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable1_ts)

# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Variable 1

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var1 <- stl(variable1_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var1 <- data.frame(
  Time = rep(time(variable1_ts), 4),  # Tiempo repetido para cada componente (son 4 componentes)
  Value = c(stl_decomp_var1$time.series[, "seasonal"], 
            stl_decomp_var1$time.series[, "trend"], 
            stl_decomp_var1$time.series[, "remainder"], 
            variable1_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable1_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var1, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 1",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

Variable 2

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var2 <- stl(variable2_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var2 <- data.frame(
  Time = rep(time(variable2_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var2$time.series[, "seasonal"], 
            stl_decomp_var2$time.series[, "trend"], 
            stl_decomp_var2$time.series[, "remainder"], 
            variable2_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable2_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var2, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 2",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

Variable 3

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var3 <- stl(variable3_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var3 <- data.frame(
  Time = rep(time(variable3_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var3$time.series[, "seasonal"], 
            stl_decomp_var3$time.series[, "trend"], 
            stl_decomp_var3$time.series[, "remainder"], 
            variable3_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable3_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 3",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

Variable 4

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var4 <- stl(variable4_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var4 <- data.frame(
  Time = rep(time(variable4_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var4$time.series[, "seasonal"], 
            stl_decomp_var4$time.series[, "trend"], 
            stl_decomp_var4$time.series[, "remainder"], 
            variable4_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable4_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var4, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 4",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

Graficamos Tasa de Crecimiento anual

Graficamos la tasa de crecimiento de la serie observada junto con la tasa de crecimiento de su tendencia. Esto nos permitió distinguir qué parte de los movimientos de cada variable corresponde a fluctuaciones transitorias de corto plazo y qué parte refleja un cambio de tendencia más estructural.

En el caso de la serie original con su tendencia, se observa con claridad que las exportaciones de pollo describen una trayectoria de crecimiento estructural con tres fases: (i) una expansión sostenida entre 2012 y 2019, (ii) el choque transitorio de 2020 seguido de una recuperación rápida, y (iii) una fase de crecimiento más moderado desde 2022, con la tendencia prácticamente plana durante 2024. Esta última fase es la que más nos interesó para efectos del pronóstico, porque sugiere que el crecimiento acelerado de la década anterior se ha agotado.

Tendencia Variable 1

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var1)

Tendencia Variable 2

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 2: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 2") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var2)

Tendencia Variable 3

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 3: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 3") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var3)

Tendencia Variable 4

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable4_vec <- as.numeric(variable4_ts)
tendencia_var4 <- as.numeric(stl_decomp_var4$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable4_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var4 <- ggplot() +
  geom_line(aes(x = fechas, y = variable4_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var4, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 4: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 4") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var4)

Tasas de crecimiento vs Tendencia

Variable 1

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var1 <- (variable1_ts[(13:length(variable1_ts))] / variable1_ts[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1 <- (tendencia_var1[(13:length(tendencia_var1))] / tendencia_var1[1:(length(tendencia_var1) - 12)] - 1) * 100

# Crear vector de fechas corregido, es decir que inicie desde enero 2013
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# Verificar longitudes
print(length(fechas_corregidas_var1))
## [1] 144
print(length(tasa_crecimiento_var1))
## [1] 144
print(length(tasa_tendencia_var1))
## [1] 144
library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable1: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

Variable 2

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))

# Verificar longitudes
print(length(fechas_corregidas_var2))
## [1] 144
print(length(tasa_crecimiento_var2))
## [1] 144
print(length(tasa_tendencia_var2))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable2: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

Variable 3

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))

# Verificar longitudes
print(length(fechas_corregidas_var3))
## [1] 144
print(length(tasa_crecimiento_var3))
## [1] 144
print(length(tasa_tendencia_var3))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var3 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable3: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

Variable 4

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var4 <- (variable4_ts[(13:length(variable4_ts))] / variable4_ts[1:(length(variable4_ts) - 12)] - 1) * 100
tasa_tendencia_var4 <- (tendencia_var4[(13:length(tendencia_var4))] / tendencia_var4[1:(length(tendencia_var4) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var4 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var4))

# Verificar longitudes
print(length(fechas_corregidas_var4))
## [1] 144
print(length(tasa_crecimiento_var4))
## [1] 144
print(length(tasa_tendencia_var4))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var4 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var4, y = tasa_crecimiento_var4), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var4, y = tasa_tendencia_var4), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable4: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var4)

Modelo ARIMA

Para este caso utilizamos el modelo ARIMA de manera manual con el fin de identificar los valores de 𝑝,𝑑,𝑞y seleccionar aquel que ofreciera la mejor capacidad predictiva o que más se ajustara al análisis que queríamos realizar. Al aplicar la prueba de Dickey-Fuller, obtuvimos un p-value menor a 0,05, lo que indica que la serie es estacionaria.

# Esta división idealmente podria se 80%-70% de los datos para entrenamiento y 20%-30% para prueba o test

# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2024 y  el conjunto de prueba o test: noviembre 2024-diciembre 2024 

train_size <- length(variable4_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable4_ts, end = c(2024, 9))  # Entrenamiento hasta septiembre 2024
test_ts <- window(variable4_ts, start = c(2024, 10))  # Prueba inicia desde oct2024

Modelo ARIMA manual

# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts)
## Warning in adf.test(train_ts): p-value smaller than printed p-value
print(adf_test)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -4.4271, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
if (adf_test$p.value > 0.05 && length(train_ts) > 1) {
  train_diff <- diff(train_ts, differences = 1)
}
# Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(FECHA = time(train_ts), POLLO = as.numeric(train_ts)), aes(x = FECHA, y = POLLO)) +
    geom_line(color = "blue") +
    ggtitle("POLLO:Serie Original") +
    xlab("FECHA") + ylab("POLLO")
  
  ggplotly(p2)  # Convertir en gráfico interactivo
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Modelo ACF

Como se observa a continuación, utilizamos para q, los valores 1,2,3,4,5,6 y para d, utilizamos 1.

# Graficar ACF y PACF
acf_plot <- ggAcf(train_ts, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_ts, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")

ggplotly(acf_plot)
ggplotly(pacf_plot)

Cálculo manual de modelo ARIMA

Como podemos observar, al realizar nuestro modelo arima con los valores(1.1.6), nos da una significancia o porcentaje funcional del modelo del 3,23, lo que quiere decir que es un modelo casi perfecto ya que está por debajo del 10%.

# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,1,6))
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(1,1,6) 
## 
## Coefficients:
##           ar1     ma1      ma2      ma3      ma4      ma5     ma6
##       -0.9207  0.8931  -0.3875  -0.3889  -0.3445  -0.1667  0.1489
## s.e.   0.1109  0.1259   0.1095   0.1073   0.1316   0.0949  0.1167
## 
## sigma^2 = 31633329:  log likelihood = -1525.31
## AIC=3066.62   AICc=3067.62   BIC=3090.81
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 1028.027 5475.335 4152.837 0.709467 3.239282 0.5282454 -0.04318113
checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,6)
## Q* = 53.582, df = 17, p-value = 1.155e-05
## 
## Model df: 7.   Total lags used: 24

Pronóstico con el modelo ARIMA manual

El primer pronóstico que realizamos se hace comparando la serie original con nuestra proyección para identificar que tan confiable es a la hora de pronosticar.

#Pronóstico con el modelo ARIMA manual
manual_forecast <- forecast(manual_arima_model, h = length(test_ts))

# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), 
                                   Pronostico = as.numeric(manual_forecast$mean),
                                   Observado = as.numeric(test_ts))

# Graficar pronóstico manual junto con los valores observados reales
p_manual <- ggplot(manual_forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico Manual")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Microempresas:Pronóstico Manual vs Observado") +
  xlab("Tiempo") + ylab("Microempresas")

ggplotly(p_manual)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# Cargar librerías necesarias
library(forecast)
library(dplyr)
## 
## 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
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_manual <- forecast(manual_arima_model, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table_manual <- data.frame(
  Tiempo = time(arima_forecast_manual$mean),  # Extraer las fechas del pronóstico
  Observado = as.numeric(test_ts),  # Valores reales
  Pronosticado = as.numeric(arima_forecast_manual$mean)  # Valores pronosticados
)

# Mostrar la tabla
print(forecast_table_manual)
##     Tiempo Observado Pronosticado
## 1 2024.750  151683.9     155241.4
## 2 2024.833  158293.7     155669.4
## 3 2024.917  156027.5     153883.6
# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo trimestre
next_month_forecast_manual <- data.frame(
  Tiempo = time(next_forecast_manual$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast_manual$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast_manual)
##     Tiempo Pronostico
## 1 2024.750   155241.4
## 2 2024.833   155669.4
## 3 2024.917   153883.6
## 4 2025.000   152337.4
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_manual, 1)
print(paste("Pronóstico para el primer mes 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para el primer mes 2025: 2025 = 152337.353606851"

SARIMA automatico

Además de nuestro modelo manual, utilizamos la función de identificación automática auto.arima(), que busca sistemáticamente la combinación de parámetros que minimiza un criterio de información estadístico. Hicimos esto para contrastar nuestra especificación manual —basada en nuestra lectura de las gráficas ACF/PACF— contra una búsqueda más sistemática y así tener un punto de comparación adicional.

# Identificación automática modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(2,1,2)(0,0,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2    sma1     drift
##       0.0776  0.2870  -0.1814  -0.7484  0.4582  410.2079
## s.e.  0.1455  0.1362   0.1098   0.1107  0.0855   77.2932
## 
## sigma^2 = 26303091:  log likelihood = -1512.98
## AIC=3039.96   AICc=3040.74   BIC=3061.13
# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(2, 1, 2),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(0, 0, 1),  # (P,D,Q) -> (1,0,0)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(2,1,2)(0,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2    sma1
##       0.0472  0.2504  -0.0957  -0.6712  0.4657
## s.e.  0.1517  0.1405   0.1166   0.1092  0.0840
## 
## sigma^2 = 27592965:  log likelihood = -1516.94
## AIC=3045.88   AICc=3046.46   BIC=3064.03
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 820.2403 5148.873 3854.356 0.5551843 3.010244 0.4902783 -0.0145149

El modelo SARIMA(2,1,2)(0,0,1)[12] ajustado muestra un buen desempeño al capturar la dinámica de la serie, ya que los coeficientes estimados son estables y el componente estacional de media móvil resulta relevante. Las métricas de ajuste, como el AIC y BIC, permiten comparar con otros modelos, mientras que los errores de entrenamiento reflejan un RMSE elevado en valores absolutos pero un MAPE cercano al 3%, lo que indica una capacidad predictiva adecuada en términos relativos.

Pronostico SARIMA AUTO

A diferencia de la gráfica anterior, se observa que la capacidad de pronóstico del modelo SARIMA es superior al ARIMA, ya que logra seguir más de cerca el comportamiento real de la serie y reduce la brecha entre los valores observados y los estimados. Esto se traduce en una mayor precisión y confiabilidad en los resultados, lo cual es especialmente útil para anticipar tendencias y tomar decisiones basadas en proyecciones más ajustadas a la realidad. Además, el componente estacional incorporado en SARIMA permite capturar patrones recurrentes que el ARIMA no logra reflejar con la misma eficacia, mejorando así la calidad del análisis y la utilidad práctica del modelo en escenarios donde la estacionalidad juega un papel importante.

# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts))  # Predecir los valores futuros

# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean), 
                            Pronostico = as.numeric(forecast_arima$mean),
                            Observado = as.numeric(test_ts))

# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("Unidad Variable 1")

ggplotly(p4)  # Convertir el gráfico en interactivo
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# Cargar librerías necesarias
library(forecast)
library(dplyr)

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(auto_arima_model, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
  Tiempo = time(arima_forecast$mean),  # Extraer las fechas del pronóstico
  Observado = as.numeric(test_ts),  # Valores reales
  Pronosticado = as.numeric(arima_forecast$mean)  # Valores pronosticados
)

# Mostrar la tabla
print(forecast_table)
##     Tiempo Observado Pronosticado
## 1 2024.750  151683.9     153599.7
## 2 2024.833  158293.7     158675.9
## 3 2024.917  156027.5     157874.6
# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
  Tiempo = time(next_forecast$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast)
##     Tiempo Pronostico
## 1 2024.750   153599.7
## 2 2024.833   158675.9
## 3 2024.917   157874.6
## 4 2025.000   154761.1
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 154761.078526839"