¿Qué es el modelo ARIMA?

Un modelo ARIMA (Autoregressive Integrated Moving Average) es una herramienta estadística utilizada para analizar y predecir series de tiempo. En términos simples, es como una “bola de cristal matemática” que usa datos pasados para estimar valores futuros, especialmente útil en finanzas para preveer precios, ventas o ingresos.

Desglose del nombre ARIMA 1. Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro 2. Integrated (I) → Ajusta tendencias en los datos para hacerlos estacionarios (sin patrones cambiantes en el tiempo). 3.MAving Average o media móvil (MA)- → Suaviza fluctuaciones aleatorias a partir de errores pasados.

ARIMA combina estos tres elementos para crear una predicción más precisa.

Metodología Box-Jenkins

La metodología Box-Jenkins es un enfoque sistemático para construir modelos ARIMA con el objetivo de analizar y pronosticar series de tiempo. Fue desarrollada por George Box y Gwilym Jenkins y se basa en cuatro etapas clave:

1️⃣ Identificación 2️⃣ Estimación 3️⃣ Validación 4️⃣ Pronóstico

Se usa especialmente cuando se quiere encontrar el modelo ARIMA más adecuado para una serie de tiempo.

Paso 1. Idntificación del Modelo

Análisis exploratorio de la variable o extracción de señales

Antes de ajustar un modelo ARIMA, es necesario entender la estructura de la serie de tiempo.

1. Verificar estacionariedad (media y varianza constantes en el tiempo) 2. Diferenciar los datos por si hay una tendencia para hacer los datos estacionarios 3. Examinar los gráficos ACF y PACF para determinar los parámetro p y q

#Cargar librerías necesarias
library(readxl)  # Para leer archivos Excel
library(tseries)  # Para pruebas de estacionariedad
library(forecast)  # Para modelado ARIMA y pronósticos
library(ggplot2)  # Para visualización de datos
library(plotly)  # Para gráficos interactivos
library(timetk)   
data_col <- read_excel("C:/Users/julie/OneDrive - PUJ Cali/Clase Predicción Econ/2025-1/Clases/Módulo III y IV/MODELO ARIMA/Colombina.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))
# Convertir los ingresos a una serie de tiempo trimestral
ingresos_ts <- ts(data_col$INGRESOS, start = c(2016, 1), frequency = 4)
ingresos_ts
##           Qtr1      Qtr2      Qtr3      Qtr4
## 2016 400871238 398002762 480807000 469500000
## 2017 397086000 398514000 452923000 478666000
## 2018 415790000 418376000 468994000 507129000
## 2019 433580000 462186000 517285000 528268000
## 2020 498207000 402864000 498107000 526615000
## 2021 485736000 462864000 580319000 632743000
## 2022 619157000 669628000 787278000 858631000
## 2023 839495000 797503000 848018000 859548000
## 2024 748846000 763244000 871052000
# Calcular estadísticas descriptivas básicas
descriptive_stats <- data.frame(
  Min = min(ingresos_ts),
  Max = max(ingresos_ts),
  Media = mean(ingresos_ts),
  Mediana = median(ingresos_ts),
  DesviacionEstandar = sd(ingresos_ts),
  CoefVar = sd(ingresos_ts) / mean(ingresos_ts)
)
print(descriptive_stats)
##         Min       Max     Media   Mediana DesviacionEstandar   CoefVar
## 1 397086000 871052000 570795229 498207000          163091152 0.2857262
# Gráfico interactivo de la serie original
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2016-01-01"), by = "quarter", length.out = nrow(data_col)), y = ingresos_ts)) +
  geom_line(color = "grey", size = 1) +
  geom_point(color = "black") +
  ggtitle("Figura 1.Ingresos Totales Trimestrales de Colombina S.A.") +
  xlab("Tiempo") +
  ylab("Miles de Pesos Colombianos") +
  theme_minimal()
ggplotly(grafico_serie)
# Extraer señales:
# Descomposición de la serie temporal
stl_decomp <- stl(ingresos_ts, s.window = "periodic")
plot(stl_decomp)

Graficamos serie original VS ajustada por estacionalidad

# Extraer los componentes de la descomposición
ingresos_sa <- ingresos_ts - stl_decomp$time.series[, "seasonal"]
library(ggplot2)
library(plotly)
library(scales)

# Comparación de la serie original con la ajustada por estacionalidad
grafico_ajustada <- ggplot() +
  geom_line(aes(x = seq_along(ingresos_ts), y = ingresos_ts), color = "grey", size = 1, linetype = "solid") +
  geom_line(aes(x = seq_along(ingresos_sa), y = ingresos_sa), color = "black", size = 1, linetype = "dashed") +
  ggtitle("Figura 2.Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Miles de Pesos Colombianos") +
  theme_minimal()
ggplotly(grafico_ajustada)
# Crear vector de fechas correctamente alineado con la serie
fechas <- seq.Date(from = as.Date("2016-01-01"), by = "quarter", length.out = length(ingresos_ts))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada <- ggplot() +
  geom_line(aes(x = fechas, y = ingresos_ts), color = "grey", size = 1, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas, y = ingresos_sa), color = "black", size = 1, linetype = "dashed", name = "Serie Ajustada") +
  ggtitle("Figura 2. Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Miles de Pesos Colombianos") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada)
# Comparación de la serie original con la tendencia
tendencia <- stl_decomp$time.series[, "trend"]
grafico_tendencia <- ggplot() +
  geom_line(aes(x = seq_along(ingresos_ts), y = ingresos_ts), color = "grey", size = 1) +
  geom_line(aes(x = seq_along(tendencia), y = tendencia), color = "purple", size = 1, linetype = "dashed") +
  ggtitle("Figura 3. Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Miles de Pesos Colombianos") +
  theme_minimal()
ggplotly(grafico_tendencia)
# Extraer la tendencia de la descomposición STL
tendencia <- stl_decomp$time.series[, "trend"]

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia <- ggplot() +
  geom_line(aes(x = fechas, y = ingresos_ts), color = "grey", size = 1, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas, y = tendencia), color = "purple", size = 1, linetype = "dashed", name = "Tendencia") +
  ggtitle("Figura 3. Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Miles de Pesos Colombianos") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_tendencia)
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento <- (ingresos_ts[(5:length(ingresos_ts))] / ingresos_ts[1:(length(ingresos_ts) - 4)] - 1) * 100
tasa_tendencia <- (tendencia[(5:length(tendencia))] / tendencia[1:(length(tendencia) - 4)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas <- seq(from = as.Date("2017-01-01"), by = "quarter", length.out = length(tasa_crecimiento))

# Verificar longitudes
print(length(fechas_corregidas))
## [1] 31
print(length(tasa_crecimiento))
## [1] 31
print(length(tasa_tendencia))
## [1] 31
# Gráfico de la tasa de crecimiento anual
grafico_crecimiento <- ggplot() +
  geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento), color = "grey", size = 1) +
  geom_line(aes(x = fechas_corregidas, y = tasa_tendencia), color = "green", size = 1, linetype = "dashed") +
  ggtitle("Figura 4. Tasa de Crecimiento Anual: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento)
# Gráfico de la tasa de crecimiento anual
grafico_crecimiento <- ggplot() +
  geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento), color = "grey", size = 1) +
  geom_line(aes(x = fechas_corregidas, y = tasa_tendencia), color = "green", size = 1, linetype = "dashed") +
  ggtitle("Figura 4. Tasa de Crecimiento Anual: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento)

Ahora dividimos el conjunto de datos de prueba y entrenamiento

✅ Para entrenar el modelo con datos históricos sin usar información futura. ✅ Para evaluar la precisión del modelo comparando sus predicciones con los datos reales de prueba.

💡 Esta división es clave en modelos predictivos para evitar sobreajuste y evaluar el rendimiento en datos no vistos.

# Dividir la serie en conjunto de entrenamiento (1T2016-3T2023) y prueba (1T2024-3T2024)
train_size <- length(ingresos_ts) - 4
train_ts <- window(ingresos_ts, end = c(2023, 3))  # Entrenamiento hasta 2023Q3
test_ts <- window(ingresos_ts, start = c(2024, 1))  # Prueba desde 2024Q1

Identificar estacionariedad

En el análisis de series de tiempo, una serie es estacionaria si su comportamiento es constante a lo largo del tiempo, es decir:

✅ Su promedio no cambia con el tiempo. ✅ Su variabilidad (qué tanto fluctúa) se mantiene estable. ✅ Su relación con valores pasados es siempre la misma.

¿Cómo saber si una serie es estacionaria?

1️⃣ Observando un gráfico Si el gráfico de la serie muestra una tendencia creciente o decreciente, o si las variaciones se hacen más grandes con el tiempo, la serie probablemente no es estacionaria.

2️⃣ Usando la prueba de Dickey-Fuller Aumentada (ADF) Es una prueba estadística que nos dice si la serie tiene una tendencia fuerte. Si el p-valor de la prueba es mayor a 0.05, significa que la serie no es estacionaria.

¿Cómo hacer una serie estacionaria?

Si encontramos que la serie no es estacionaria, podemos transformarla para que lo sea:

✅ Diferenciación: Restamos cada valor con su valor anterior. Esto elimina tendencias crecientes o decrecientes. ✅ Tomar el logaritmo: Si la variabilidad crece con el tiempo, aplicar un logaritmo estabiliza la varianza. ✅ Eliminar tendencias o ajustar estacionalidad: Si hay patrones repetitivos, podemos restarlos o modelarlos por separado.

# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts)
print(adf_test)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -0.068677, Lag order = 3, p-value = 0.99
## alternative hypothesis: stationary
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
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(Tiempo = time(train_ts), Ingresos = as.numeric(train_ts)), aes(x = Tiempo, y = Ingresos)) +
    geom_line(color = "blue") +
    ggtitle("Serie Original") +
    xlab("Tiempo") + ylab("Ingresos")
  
  ggplotly(p2)  # Convertir en gráfico interactivo
  # Graficar la serie diferenciada en un gráfico separado
  p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], Ingresos_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = Ingresos_Diff)) +
    geom_line(color = "red") +
    ggtitle("Serie Estacionaria (Diferenciada)") +
    xlab("Tiempo") + ylab("Ingresos Diferenciados")
  
  ggplotly(p3)  # Convertir en gráfico interactivo
# Segunda prueba de estacionariedad sobre la serie diferenciada
  adf_test_diff <- adf.test(train_diff)
  print(adf_test_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff
## Dickey-Fuller = -1.8608, Lag order = 3, p-value = 0.6262
## alternative hypothesis: stationary

A continuación se muestra que el valor que puede tomar d=2:

El valor p ya es menor a 0.05.

train_diff2 <- diff(train_ts, differences = 2)  # Segunda diferenciación
adf_test_diff2 <- adf.test(train_diff2)
print(adf_test_diff2)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff2
## Dickey-Fuller = -4.6847, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Identificación manual de p y q

¿Qué hacen estos gráficos?

ACF (Autocorrelation Function)

Muestra la correlación de la serie con sus rezagos.

Ayuda a determinar el parámetro q en un modelo ARIMA(p, d, q).

PACF (Partial Autocorrelation Function)

Muestra la correlación parcial entre la serie y un rezago específico, eliminando el efecto de rezagos intermedios.

Ayuda a determinar el parámetro p en un modelo ARIMA(p, d, q).

# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff2, lag.max = 20) + ggtitle("Autocorrelation Function (ACF)")
pacf_plot <- ggPacf(train_diff2, lag.max = 20) + ggtitle("Partial Autocorrelation Function (PACF)")

ggplotly(acf_plot)
ggplotly(pacf_plot)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(0,2,1))
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.9999
## s.e.   0.2476
## 
## sigma^2 = 3.269e+15:  log likelihood = -560.33
## AIC=1124.65   AICc=1125.11   BIC=1127.39
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 6163197 54337061 44717139 0.3434413 8.543093 0.6550887 0.02872554

Validación de residuos

 checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,1)
## Q* = 39.946, df = 5, p-value = 1.531e-07
## 
## Model df: 1.   Total lags used: 6
#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("Pronóstico Manual vs Observado") +
  xlab("Tiempo") + ylab("Ingresos")

ggplotly(p_manual)
# Calcular métricas de evaluación del modelo manual
mae_manual <- mean(abs(manual_forecast$mean - test_ts), na.rm = TRUE)
rmse_manual <- sqrt(mean((manual_forecast$mean - test_ts)^2, na.rm = TRUE))

# Mostrar métricas de evaluación del modelo manual
cat("MAE Manual: ", mae_manual, "\n")
## MAE Manual:  129235352
cat("RMSE Manual: ", rmse_manual, "\n")
## RMSE Manual:  129235601
# Identificación automática del mejor modelo ARIMA
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(0,1,0)(0,1,1)[4] 
## 
## Coefficients:
##          sma1
##       -0.7176
## s.e.   0.2572
## 
## sigma^2 = 1.355e+15:  log likelihood = -490.76
## AIC=985.53   AICc=986.05   BIC=988.04

ARIMA(0,1,0)(0,1,1)[4]

🔹 Desglose de los parámetros:

ARIMA(0,1,0) → No tiene términos autorregresivos (AR), tiene una diferenciación (d=1), y no tiene términos de media móvil (MA). (0,1,1)[4] → Es un modelo estacional, con: 0 términos autorregresivos estacionales. 1 diferenciación estacional. 1 término de media móvil estacional. Periodicidad 4, lo que sugiere un patrón trimestral.

El modelo anterior sugiere que los ingresos tienen un fuerte componente estacional trimestral, y se necesita una diferenciación regular y estacional para hacer la serie estacionaria.

# Ajuste del modelo ARIMA con los parámetros encontrados
darima <- Arima(train_ts, order = auto_arima_model$arma[1:3])  # Ajustar modelo con los parámetros identificados
summary(darima)
## Series: train_ts 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##            mean
##       539843323
## s.e.   26009382
## 
## sigma^2 = 2.108e+16:  log likelihood = -626.08
## AIC=1256.16   AICc=1256.59   BIC=1259.03
## 
## Training set error measures:
##                        ME      RMSE       MAE       MPE     MAPE     MASE
## Training set 7.688575e-08 142827109 114463361 -5.919822 20.54873 1.676844
##                   ACF1
## Training set 0.8320102
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with non-zero mean
## Q* = 64.918, df = 6, p-value = 4.484e-12
## 
## Model df: 0.   Total lags used: 6
# 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("Ingresos")

ggplotly(p4)  # Convertir el gráfico en interactivo