Series de tiempo: Extracción de señales y Modelo ARIMA
Julieth Cerón
2025-03-12
Introducción
En este informe, se realizará un análisis de series temporales: *1) Extracción de señales o descomposición de series temporales y 2) Pronóstico del número de microempresas nuevas en Cali. El periodo de análisis es desde el primer trimestre de 2004 (1T2004) y el cuarto trimestre de 2024 (4T2024).
El análisis de series temporales permite identificar patrones subyacentes en los datos, como tendencias, estacionalidad y componentes irregulares, lo que facilita una mejor comprensión del comportamiento de la variable a analizar. Para ello, se emplearán técnicas estadísticas o métodos de descomposición y el Modelo ARIMA con el fin de extraer información clave que pueda contribuir a la toma de decisiones estratégicas.
A lo largo del informe, se presentarán los resultados obtenidos y se discutirán sus implicaciones en el contexto empresarial de Cali.
Puntualmente, el análisis de la creación de microempresas en Cali es clave para la planificación económica y la toma de decisiones estratégicas. Identificar patrones y tendencias en estos datos ayuda a mejorar las condiciones para los emprendedores y fortalecer el desarrollo económico local.
Análisis exploratorio de la variable
#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) Cargar base de datos
library(readxl)
data_col <- read_excel("MICRO.xlsx", col_types = c("date",
"numeric"))Un paso indispensable es declarar la variable como serie temporal:
# Convertir/declarar el número de microempresas en serie de tiempo trimestral
micro_ts <- ts(data_col$MICRO, start = c(2004, 1), frequency = 4)
micro_ts## Qtr1 Qtr2 Qtr3 Qtr4
## 2004 4064 3394 3282 1951
## 2005 2956 3188 3768 1812
## 2006 3436 3241 3267 2058
## 2007 4327 3408 3255 2168
## 2008 3592 3300 3489 2090
## 2009 3281 3027 3328 2340
## 2010 3503 3205 3469 2006
## 2011 3493 3503 3457 2131
## 2012 3847 3412 3447 2224
## 2013 3462 3739 3896 2715
## 2014 4406 3686 3900 3122
## 2015 3609 3460 3784 2459
## 2016 4299 4072 4276 2769
## 2017 5106 3896 3956 2618
## 2018 5009 4406 4349 3481
## 2019 6508 5169 5488 3510
## 2020 6046 2600 6156 4396
## 2021 6197 3860 4803 3956
## 2022 6210 4869 5785 3700
## 2023 5362 4409 4485 3087
## 2024 4939 4530 4473 3150Ahora si calculamos estadísticas descriptivas
En la tabla siguiente se observa lo siguiente:
En promedio, cada trimestre se han creado alrededor de 3760 microempresas. Este es un valor de referencia para entender el comportamiento típico.
La mediana indica que la mitad de los trimestres tuvieron menos de 3506.5 microempresas nuevas y la otra mitad tuvo más. Como está un poco por debajo de la media, sugiere que hay algunos valores altos que están elevando el promedio.
Existe una variabilidad notable en la cantidad de microempresas nuevas por trimestre. Un valor de 1052 indica que los datos fluctúan bastante alrededor del promedio.
la variabilidad relativa de la serie no es excesivamente alta. Un 27.99% implica que, aunque hay oscilaciones, los datos no son extremadamente dispersos respecto al promedio.
# Calcular estadísticas descriptivas básicas
descriptive_stats <- data.frame(
Min = min(micro_ts),
Max = max(micro_ts),
Media = mean(micro_ts),
Mediana = median(micro_ts),
DesviacionEstandar = sd(micro_ts),
CoefVar = sd(micro_ts) / mean(micro_ts)
)
print(descriptive_stats)## Min Max Media Mediana DesviacionEstandar CoefVar
## 1 1812 6508 3760.5 3506.5 1052.423 0.2798625Exploramos el comportamiento o evolución de la variable:
En la Figura 1. destaca lo siguiente:
- A partir de 2015 el número de microempresas nuevas aumenta significativamente. Sin embargo, a medida que el número de microempresas aumenta, también crecen las fluctuaciones entre trimestres. Hay períodos con picos muy altos y caídas pronunciadas.
*Hacia el final de la serie se nota una caída en el número de microempresas. Podría estar asociado a factores económicos como inflación, desaceleración económica o menor acceso a financiamiento para emprendedores.
Por tanto, la Figura 1, sugiere que la creación de microempresas en Cali ha aumentado a lo largo del tiempo, pero con fuertes fluctuaciones. Es clave analizar qué factores han impulsado los picos y qué ha causado las caídas recientes, para entender mejor el ecosistema emprendedor de la ciudad.
# Gráfico interactivo de la serie original
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2004-01-01"), by = "quarter", length.out = nrow(data_col)), y = micro_ts)) +
geom_line(color = "grey", size = 0.4) +
geom_point(color = "black", size = 0.1) +
ggtitle("Figura 1.Evolución de las microempresas en Cali") +
xlab("Tiempo") +
ylab("Número de microempresas") +
theme_minimal()
ggplotly(grafico_serie)Extracción de señales
*Muchas series de tiempo son una combinación de varias influencias. Es por eso que, separar la tendencia, la estacionalidad y los componentes aleatorios permite entender mejor qué está impulsando los cambios en la serie.
Si analizamos la creación de microempresas en Cali, podríamos querer saber si el crecimiento se debe a una tendencia real o a fluctuaciones estacionales.
- Los modelos de pronóstico funcionan mejor cuando las señales subyacentes están bien definidas. Por ejemplo, si eliminamos la estacionalidad de una serie financiera, los modelos predictivos pueden enfocarse en la tendencia real y reducir errores.
*Detectar cambios inesperados en la serie es más fácil cuando se eliminan componentes predecibles. Ejemplo: Si hay una caída abrupta en el número de microempresas, podemos verificar si es una anomalía (ruido) o un cambio estructural en la economía.
En conclusión, la descomposición de series de tiempo permite comprender mejor los datos, mejorar predicciones y tomar decisiones más estratégicas. Es una herramienta clave en la analítica de negocios, especialmente en entornos donde las fluctuaciones en los datos pueden afectar inversiones, políticas económicas y estrategias empresariales.
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp <- stl(micro_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df <- data.frame(
Time = rep(time(micro_ts), 4), # Tiempo repetido para cada componente
Value = c(stl_decomp$time.series[, "seasonal"],
stl_decomp$time.series[, "trend"],
stl_decomp$time.series[, "remainder"],
micro_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(micro_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Figura A. Descomposición del número de microempresas en Cali",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Interpretación Figura A.Descomposición temporal
Componente estacional
Muestra patrones recurrentes a lo largo del tiempo. La interpretación de este componente es clave para entender fluctuaciones predecibles dentro de cada período analizado.
Las caidas o incrementos del gcomponente estacional ocurren a intervalos similares (Primer trimestre y cuarto trimestre de cada año), lo que indica que el comportamiento estacional se mantiene a lo largo del tiempo.
Así, el componente estacional sugiere que hay fluctuaciones repetitivas en la creación de microempresas en Cali a lo largo de cada año. Comprender estos patrones permite ajustar estrategias de apoyo al emprendimiento y realizar pronósticos más precisos para futuras decisiones económicas. También puede ayudar a prever en qué períodos del año se espera un menor número de empresas y planificar estrategias para mitigar estos efectos.
Graficamos serie original VS ajustada por estacionalidad
En la Figura 2, se observa que en la serie ajustada por estacionalidad se pueden apreciar con mayor precisión cambios estructurales, facilitando el análisis sin la interferencia de variaciones temporales recurrentes.
Estos análisis permiten diseñar estrategias de apoyo al emprendimiento basadas en datos reales, sin sesgos estacionales. Además, proporciona información clave para la toma de decisiones estratégicas en el sector de microempresas y, mejora la precisión de las estimaciones futuras al eliminar patrones estacionales que podrían distorsionar los análisis.
En conclusión, la descomposición de la serie temporal no solo facilita la identificación de tendencias de largo plazo, sino que también permite tomar decisiones más informadas, anticipando oportunidades y riesgos en el sector de las microempresas en Cali
# Extraer los componentes de la descomposición
micro_sa <- micro_ts - stl_decomp$time.series[, "seasonal"]# Crear vector de fechas correctamente alineado con la serie
fechas <- seq.Date(from = as.Date("2004-01-01"), by = "quarter", length.out = length(micro_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada <- ggplot() +
geom_line(aes(x = fechas, y = micro_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas, y = micro_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Figura 2. Microempresas:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Número de microempresas") +
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)Ahora graficamos serie original vs tendencia
La serie original de microempresas contiene variaciones estacionales y aleatorias que pueden dificultar la identificación de patrones reales. La extracción de la tendencia permite centrarse en los cambios estructurales de la serie.
Analizar la tendencia ayuda a prever escenarios futuros y anticipar posibles crisis o oportunidades en el sector de microempresas.
*Después del fuerte crecimiento post-2020, la tendencia parece estabilizarse e incluso mostrar una leve reducción, sugiriendo posibles cambios en las condiciones económicas o en el mercado.
# 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 = micro_ts), color = "grey", size = 0.7, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas, y = tendencia), color = "black", size = 0.8, linetype = "solid", name = "Tendencia") +
ggtitle("Figura 3. Microempresas:Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Número de microempresas") +
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)Ahora calculamos la tasa de crecimiento de la serie original vs tendencia
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento <- (micro_ts[(5:length(micro_ts))] / micro_ts[1:(length(micro_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("2005-01-01"), by = "quarter", length.out = length(tasa_crecimiento))
# Verificar longitudes
print(length(fechas_corregidas))## [1] 80print(length(tasa_crecimiento))## [1] 80print(length(tasa_tendencia))## [1] 80# Gráfico de la tasa de crecimiento anual
grafico_crecimiento <- ggplot() +
geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas, y = tasa_tendencia), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Figura 4. Microempresas: Tasa de crecimiento anual % de la serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento)Interpretación Figura 3. Tendencia VS original: crecimiento anual%
Mientras que la serie original muestra mucha variabilidad, la tendencia permite identificar patrones más claros y estructurales en el crecimiento de las microempresas.
Analizar la tasa de crecimiento anual ayuda a detectar cambios en el entorno económico que afectan el sector. Se pueden prever crisis o períodos de auge y prepararse para ellos.
Permite evaluar si el sector de microempresas en Cali sigue una dinámica similar a la nacional o global.
*A lo largo del tiempo, la tasa de crecimiento de las microempresas ha sido volátil, alternando entre períodos de expansión y contracción.
*Se identifican picos de crecimiento antes de 2020 y una fuerte caída alrededor de ese año, probablemente por la pandemia.
*La caída abrupta en 2020 es seguida por una recuperación rápida, lo que indica la resiliencia del sector.
*Sin embargo, la tendencia reciente sugiere menor estabilidad en el crecimiento de nuevas microempresas.
¿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
- Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro
- 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.
Antes de empezar a aplicar la metodlogía BOX-JENKINS, debemos dividir 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 (1T2004-1T2024) y prueba (2T2024-4T2024)
train_size <- length(micro_ts) - 2
train_ts <- window(micro_ts, end = c(2024, 1)) # Entrenamiento hasta 2024Q1
test_ts <- window(micro_ts, start = c(2024, 2)) # Prueba desde 2024Q2Paso 1: Identificación del modelo
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.
Test de Dickey-Fuller
El test de Dickey-Fuller aumentado (ADF) se usa para verificar si una serie temporal es estacionaria, es decir, si sus propiedades estadísticas (media y varianza) permanecen constantes en el tiempo.
HO: Serie no estacionaria HI: Serie estacionaria
¿Qué significa el p-valor?
Si el p-valor es bajo (< 0.05) → Rechazamos la hipótesis nula y concluimos que la serie es estacionaria. Si el p-valor es alto (> 0.05) → No podemos rechazar la hipótesis nula, lo que indica que la serie no es estacionaria.
# 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 = -2.4893, Lag order = 4, p-value = 0.3754
## 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)
}Diferenciación en niveles
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), micro = as.numeric(train_ts)), aes(x = Tiempo, y = micro)) +
geom_line(color = "blue") +
ggtitle("Microempresas:Serie Original") +
xlab("Tiempo") + ylab("Microempresas")
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], micro_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = micro_Diff)) +
geom_line(color = "red") +
ggtitle("Microempresas:Serie Estacionaria (Una diferenciación)") +
xlab("Tiempo") + ylab("Ingresos Diferenciados")
ggplotly(p3) # Convertir en gráfico interactivo*Diferenciación en logaritmo**
Cuando una serie de tiempo tiene una creciente varianza (lo que significa que la amplitud de las fluctuaciones aumenta con el tiempo), aplicar un logaritmo puede ayudar a estabilizar esta varianza. Muchas series económicas o financieras, como el precio de acciones o el Producto Interno Bruto (PIB), tienden a mostrar crecimiento exponencial o crecimiento en porcentaje (por ejemplo, tasas de crecimiento de doble dígito).
En conclusión, la aplicación de logaritmos en series de tiempo se realiza principalmente para lograr que la serie sea más estable, lineal y estacionaria. Esta transformación es relevante porque permite modelar mejor las series que siguen un crecimiento exponencial y facilita la aplicación de técnicas estadísticas que requieren estacionariedad.
# 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 = -2.4893, Lag order = 4, p-value = 0.3754
## 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_log <- diff(log(train_ts), differences = 1)
}
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), micro = as.numeric(train_ts)), aes(x = Tiempo, y = micro)) +
geom_line(color = "blue") +
ggtitle("Microempresas:Serie Original") +
xlab("Tiempo") + ylab("Microempresas")
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], micro_Diff = as.numeric(train_diff_log)), aes(x = Tiempo, y = micro_Diff)) +
geom_line(color = "red") +
ggtitle("Microempresas:Serie Estacionaria (Una diferenciación en logaritmo)") +
xlab("Tiempo") + ylab("microempresas Diferenciados")
ggplotly(p3) # Convertir en gráfico interactivoAhora probamos estacionariedad en la serie diferenciada (nivel y logaritmo)
# 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 = -5.7016, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
adf_test_diff_log <- adf.test(train_diff_log)
print(adf_test_diff_log)##
## Augmented Dickey-Fuller Test
##
## data: train_diff_log
## Dickey-Fuller = -5.8046, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryEn el test ADF se muestra que el valor que puede tomar d=1:
El valor-p ya es menor a 0.05 con una primera diferencia en niveles o con logaritmo natural.
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).
Recordemos El eje X representa los rezagos (lags). El eje Y muestra la autocorrelación parcial en cada rezago. Las líneas azules punteadas son los intervalos de confianza (aproximadamente 95%). Si una barra sobrepasa estos límites, indica una autocorrelación significativa. Si las barras caen dentro de los límites, no son significativamente diferentes de cero
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)ggplotly(pacf_plot)Paso 2. Estimación del modelo
El modelo manual siguiente refleja:
Coeficientes AR (ar1, ar2, ar3): Todos son negativos y cercanos a -1, lo que indica una fuerte dependencia de los valores pasados.
AIC (1252.02) y BIC (1271.08): Se usan para comparar modelos; cuanto más bajos, mejor.
Métricas de evaluación
Mean Absolute Error (MAE) = 374.491
Representa el error absoluto promedio entre las predicciones del modelo y los valores reales. Interpretación: En promedio, el modelo se equivoca en 374 microempresas al predecir la cantidad de nuevas microempresas en Cali.
Root Mean Squared Error (RMSE) = 518.0126
Similar al MAE, pero da más peso a los errores grandes, porque eleva las diferencias al cuadrado antes de promediarlas. Interpretación: Un error típico en la predicción es de aproximadamente 518 microempresas.
Comparación con MAE: Como el RMSE es mayor que el MAE, es posible que haya algunos errores grandes que estén influyendo más en el RMSE.
Mean Absolute Percentage Error (MAPE) = 10.17%
Expresa el error en términos relativos, como porcentaje del valor real.
Interpretación: En promedio, el modelo se equivoca en un 10.17% al predecir el número de microempresas nuevas.
Regla general: MAPE < 10% → Muy buen modelo ✅ 10%-20% → Modelo aceptable 👍 20%-50% → Modelo pobre ⚠️ 50% → Modelo muy malo ❌
En este caso, un MAPE de 10.17% sugiere un modelo aceptable, pero no perfecto.
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(3,1,4))
summary(manual_arima_model)## Series: train_ts
## ARIMA(3,1,4)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 ma4
## -1.0164 -1.0167 -0.9958 0.1952 0.3502 0.2239 -0.5915
## s.e. 0.0130 0.0130 0.0064 0.1064 0.1132 0.1107 0.1163
##
## sigma^2 = 297744: log likelihood = -618.01
## AIC=1252.02 AICc=1254.05 BIC=1271.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 30.93174 518.0126 374.491 -0.8938534 10.17033 0.8049523 0.1075101Significancia de coefcientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.0164122 0.0129769 -78.3247 < 2.2e-16 ***
## ar2 -1.0166948 0.0129681 -78.3996 < 2.2e-16 ***
## ar3 -0.9957855 0.0063695 -156.3365 < 2.2e-16 ***
## ma1 0.1952436 0.1063745 1.8354 0.066441 .
## ma2 0.3502133 0.1131512 3.0951 0.001968 **
## ma3 0.2239080 0.1106627 2.0233 0.043038 *
## ma4 -0.5914626 0.1163049 -5.0854 3.668e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ar1, ar2 y ar3 → Son altamente significativos (***), lo que significa que estos coeficientes tienen un impacto importante en el modelo.
** Paso 3. Validación de residuos del modelo estimado manual**
La validación de residuos es crucial para determinar si el modelo ARIMA es adecuado o si necesita mejoras. El objetivo es verificar que los residuos (errores de predicción) se comporten como ruido blanco, es decir, sin patrones detectables.
1. Serie de residuos (gráfico superior)
Muestra cómo se comportan los errores a lo largo del tiempo. Idealmente, deberían oscilar alrededor de cero sin tendencias evidentes ni grandes acumulaciones de error. Problema posible: Se observan picos alrededor de 2020, lo que podría indicar cambios estructurales en los datos (como efectos de la pandemia en la creación de microempresas).
Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)
Si el modelo es adecuado, los residuos no deben mostrar correlaciones significativas en el tiempo.
Interpretación: La mayoría de las barras están dentro de las líneas azules (intervalos de confianza). Sin embargo, algunos rezagos parecen salir del rango, lo que sugiere que aún puede haber estructura no capturada en los datos. Esto indica que el modelo podría mejorarse pero es un modelo aceptable.
Histograma de residuos con ajuste normal (gráfico inferior derecho)
Sirve para verificar si los errores siguen una distribución normal, lo cual es un supuesto clave en ARIMA. Interpretación: La curva roja representa la distribución normal teórica. Los residuos se acercan a la normalidad, pero hay algunos valores extremos (colas más gruesas de lo esperado). Esto indica que puede haber eventos atípicos o datos no bien explicados por el modelo.
Conclusión y acciones recomendadas
✅ El modelo ARIMA(3,1,4) parece razonablemente bueno, pero tiene algunas señales de que podría mejorarse.
⚠️ Posibles mejoras:
Revisar la estructura del modelo: Se pueden probar otros órdenes ARIMA o incluso modelos más complejos como SARIMA o modelos con variables exógenas (ARIMAX).
Manejo de valores atípicos: Considerar incluir un término de intervención si eventos como la pandemia afectaron la serie.
Transformación de datos: Si los residuos no son normales, una transformación logarítmica o Box-Cox puede ayudar.
Recordar: El supuesto de normalidad significa que los errores o residuos de un modelo deben seguir una distribución normal (o “campana de Gauss”). Si los errores son normales, podemos hacer predicciones más confiables y usar ciertas pruebas estadísticas que asumen esta propiedad.
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,1,4)
## Q* = 7.8333, df = 3, p-value = 0.04959
##
## Model df: 7. Total lags used: 10Paso 4. Pronóstico (modelo manual)
El modelo sigue la tendencia general, pero tiene un sesgo de sobreestimación.
Si la serie tiene patrones estacionales fuertes y no están bien capturados, el modelo puede fallar en prever las fluctuaciones.
#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)# 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: 278.2542cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 348.61MAE (Mean Absolute Error)
Indica el error promedio en unidades de la variable, es decir, en número de microempresas. En promedio, el modelo se equivoca en 278 microempresas al hacer su predicción.
RMSE (Root Mean Squared Error)
Penaliza más los errores grandes debido a la elevación al cuadrado antes de calcular la raíz. Su valor de 348.61 microempresas sugiere que hay algunos errores grandes que afectan la precisión del modelo.
Tabla de pronóstico modelo manual
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# 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.25 4530 3955.905
## 2 2024.50 4473 4626.007
## 3 2024.75 3150 3257.661Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueba una observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.
# 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_quarter_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_quarter_forecast_manual)## Tiempo Pronostico
## 1 2024.25 3955.905
## 2 2024.50 4626.007
## 3 2024.75 3257.661
## 4 2025.00 5268.218# Extraer solo el valor del trimestre adicional (último de la tabla)
next_quarter <- tail(next_quarter_forecast_manual, 1)
print(paste("Pronóstico para el primer trimestre 2025:", next_quarter$Tiempo, "=", next_quarter$Pronostico))## [1] "Pronóstico para el primer trimestre 2025: 2025 = 5268.21779519335"Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual
Si no hay test y solo quieres el siguiente punto u observación, el código a usar seria algo asi:
# Pronosticar hasta el segundo trimestre de 2024
future_forecast_manual <- forecast(manual_arima_model, h = 1)
# Extraer el valor específico del el segundo trimestre de 2024
forecast_2024q2 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para el segundo trimestre de 2024:", forecast_2024q2))## [1] "Pronóstico para el segundo trimestre de 2024: 3955.90518357422"Modelo ARIMA automático
Usa la función auto.arima() de forecast en R para seleccionar automáticamente los mejores parámetros (p,d,q).
✅ Ventajas:
✔ Optimización automática: Encuentra los valores óptimos de ARIMA sin intervención manual. ✔ Ahorra tiempo: Útil cuando hay muchas series a modelar. ✔ Evita sesgo humano: Reduce el riesgo de elegir un modelo incorrecto por falta de experiencia. ✔ Incluye corrección por estacionalidad si se usa con seasonal = TRUE. ✔ Suele funcionar bien en la mayoría de los casos, ya que usa criterios como AIC/BIC para optimizar.
❌ Desventajas: ❌ Puede no ser el mejor modelo posible, ya que depende del criterio de selección. ❌ Menos interpretabilidad: No siempre es claro por qué eligió ciertos parámetros. ❌ Puede ignorar conocimiento experto sobre la serie o factores externos.
Identificación automática del modelo ARIMA
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)
# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)## Series: train_ts
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## -0.4384 -0.7675
## s.e. 0.1086 0.0698
##
## sigma^2 = 678663: log likelihood = -650.46
## AIC=1306.92 AICc=1307.23 BIC=1314.06
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 70.95191 808.4102 668.5027 -3.050944 19.81266 1.436918 0.02937838Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.438393 0.108621 -4.036 5.438e-05 ***
## ma1 -0.767478 0.069781 -10.998 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(1,1,1) sin parte estacional
darima_auto <- Arima(train_ts,
order = c(1, 1, 1)) # Especificamos directamente (p=1, d=1, q=1)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## -0.4384 -0.7675
## s.e. 0.1086 0.0698
##
## sigma^2 = 678663: log likelihood = -650.46
## AIC=1306.92 AICc=1307.23 BIC=1314.06
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 70.95191 808.4102 668.5027 -3.050944 19.81266 1.436918 0.02937838# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 49.64, df = 6, p-value = 5.551e-09
##
## Model df: 2. Total lags used: 8Pornóstico modelo ARIMA automático
# Generar pronóstico para el conjunto de prueba
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data_auto <- data.frame(Tiempo = time(forecast_arima_auto$mean),
Pronostico = as.numeric(forecast_arima_auto$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4auto <- ggplot(forecast_data_auto, 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("Microempresas")
ggplotly(p4auto) # Convertir el gráfico en interactivo# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_auto <- data.frame(
Tiempo = time(arima_forecast_auto$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_auto$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_auto)## Tiempo Observado Pronosticado
## 1 2024.25 4530 4231.445
## 2 2024.50 4473 4541.632
## 3 2024.75 3150 4405.648Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueb auna observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo trimestre
next_quarter_forecast_auto <- data.frame(
Tiempo = time(next_forecast_auto$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_auto$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_quarter_forecast_auto)## Tiempo Pronostico
## 1 2024.25 4231.445
## 2 2024.50 4541.632
## 3 2024.75 4405.648
## 4 2025.00 4465.263# Extraer solo el valor del trimestre adicional (último de la tabla)
next_quarter <- tail(next_quarter_forecast_auto, 1)
print(paste("Pronóstico para el primer trimestre 2025:", next_quarter$Tiempo, "=", next_quarter$Pronostico))## [1] "Pronóstico para el primer trimestre 2025: 2025 = 4465.26254114641"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el 2T2024
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 1)
# Extraer el valor específico del 2T2024
forecast_2024Q2_auto <- future_forecast_auto$mean[1]
print(paste("Pronóstico para el segundo trimestre de 2024:", forecast_2024Q2_auto))## [1] "Pronóstico para el segundo trimestre de 2024: 4231.44464709905"Identificación automática del modelo SARIMA
Un SARIMA (Seasonal ARIMA) es un modelo que combina:
✅ Un modelo ARIMA tradicional para capturar relaciones entre valores pasados y errores. ✅ Un componente estacional, útil cuando los datos muestran patrones repetitivos en el tiempo (por ejemplo, ventas trimestrales, datos mensuales de temperatura, etc.).
Básicamente, es como un ARIMA mejorado que puede manejar ciclos estacionales.
En este caso, el modelo detectó un patrón repetitivo cada 4 períodos, por lo que usa componentes autoregresivos y de media móvil para modelarlo
Este modelo se descompone en dos partes:
Parte No Estacional (ARIMA(p,d,q)
AR(1): Hay una relación de la serie con su primer rezago, es decir, el valor actual depende del valor anterior con cierto peso (coeficiente ar1 = 0.2492).
d = 0: No hay diferenciación en la serie, porque lo hace en diferencia estacional.
MA(0): No hay un componente de media móvil en la parte no estacional.
Parte Estacional (SARIMA(P,D,Q)[m])
SAR(1): La serie también tiene una dependencia con su valor en la misma temporada anterior (coeficiente sar1 = 0.3745).
D = 1: Se aplica una diferenciación estacional para eliminar tendencias estacionales.
SMA(1): Hay un componente de media móvil estacional, es decir, un efecto de ruido de la temporada anterior influye en la temporada actual (coeficiente sma1 = -0.9399).
m = 4: La periodicidad estacional es de 4 períodos, lo que significa que hay un patrón repetitivo cada 4 observaciones.
Drift (deriva)
drift = 24.43: Indica una tendencia lineal en la serie en el tiempo, lo que sugiere un crecimiento constante.
# 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(1,0,0)(1,1,1)[4] with drift
##
## Coefficients:
## ar1 sar1 sma1 drift
## 0.2492 0.3745 -0.9399 24.4306
## s.e. 0.1178 0.2032 0.2441 5.2394
##
## sigma^2 = 289446: log likelihood = -594.12
## AIC=1198.24 AICc=1199.08 BIC=1209.96# Ajuste del modelo ARIMA con los parámetros encontrados
darima <- Arima(train_ts,
order = c(auto_arima_model$arma[1], # p
auto_arima_model$arma[6], # d
auto_arima_model$arma[2]), # q
seasonal = list(order = c(auto_arima_model$arma[3], # P
auto_arima_model$arma[7], # D
auto_arima_model$arma[4]), # Q
period = auto_arima_model$arma[5])) # m (periodicidad estacional)
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(1,0,0)(1,1,1)[4]
##
## Coefficients:
## ar1 sar1 sma1
## 0.3299 0.1341 -0.5968
## s.e. 0.1219 0.2420 0.2061
##
## sigma^2 = 326102: log likelihood = -597.12
## AIC=1202.24 AICc=1202.8 BIC=1211.61
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 93.14452 545.821 368.9106 1.118311 9.628784 0.7929576 -0.07045355# 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(1,0,0)(1,1,1)[4]
## Q* = 3.2224, df = 5, p-value = 0.6657
##
## Model df: 3. Total lags used: 8# 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# 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.25 4530 4478.542
## 2 2024.50 4473 4884.316
## 3 2024.75 3150 3536.506Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueba una observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo trimestre
next_quarter_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_quarter_forecast)## Tiempo Pronostico
## 1 2024.25 4478.542
## 2 2024.50 4884.316
## 3 2024.75 3536.506
## 4 2025.00 5378.273# Extraer solo el valor del trimestre adicional (último de la tabla)
next_quarter <- tail(next_quarter_forecast, 1)
print(paste("Pronóstico para el primer trimestre 2025:", next_quarter$Tiempo, "=", next_quarter$Pronostico))## [1] "Pronóstico para el primer trimestre 2025: 2025 = 5378.27272900087"Otra forma para calcular un valor futuro (fuera de muestra)
# 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(1,0,0)(1,1,1)[4] with drift
##
## Coefficients:
## ar1 sar1 sma1 drift
## 0.2492 0.3745 -0.9399 24.4306
## s.e. 0.1178 0.2032 0.2441 5.2394
##
## sigma^2 = 289446: log likelihood = -594.12
## AIC=1198.24 AICc=1199.08 BIC=1209.96# Pronosticar hasta el 2T2024
future_forecast_sarima <- forecast(auto_arima_model, h = 1)
# Extraer el valor específico del primer trimestre de 2025
forecast_2024Q2_sarima <- future_forecast_sarima$mean[1]
print(paste("Pronóstico para el segundo trimestre de 2024:", forecast_2024Q2_sarima))## [1] "Pronóstico para el segundo trimestre de 2024: 4478.54178979002"