bancolombia
LUISA ERAZO
2025-03-19
INTRODUCCIÓN
En este informe, se realizará la extracción de señales o descomposición de series temporales aplicadas al precio de la acción de Bancolombia mensual, durante el período comprendido entre el 2010 y 2025
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 del precio de la acción de la empresa. Para ello, se emplearán técnicas estadísticas y métodos de descomposición, 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 financiero de Bancolombia
Series de tiempo: Extracción de señales y Modelo ARIMA
#Cargar librerías necesarias
library(readxl) # Para leer archivos Excel## Warning: package 'readxl' was built under R version 4.4.3library(tseries) # Para pruebas de estacionariedad## Warning: package 'tseries' was built under R version 4.4.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast) # Para modelado ARIMA y pronósticos## Warning: package 'forecast' was built under R version 4.4.3library(ggplot2) # Para visualización de datos## Warning: package 'ggplot2' was built under R version 4.4.3library(plotly) # Para gráficos interactivos## Warning: package 'plotly' was built under R version 4.4.3##
## 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':
##
## layoutlibrary(timetk) ## Warning: package 'timetk' was built under R version 4.4.3library(readxl)
data_col <- read_excel("C:/Users/luisi/Downloads/PACCIONES COL.xlsx", col_types = c("date", "numeric"))## New names:
## • `` -> `...1`# Convertir/declarar variable 1=BIC en serie de tiempo mensual
variable1_ts <- ts(data_col$BIC, start = c(2010, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2010 22060.15 22159.52 22139.65 22874.99 23351.96 23610.32 26134.33 28240.97
## 2011 26849.79 26730.55 29016.06 28618.58 28916.69 29016.06 28817.32 28618.58
## 2012 27980.00 27900.00 28100.00 28900.00 27200.00 26980.00 27180.00 25800.00
## 2013 30480.00 28800.00 28800.00 30000.00 27500.00 26400.00 26800.00 26380.00
## 2014 22000.00 24100.00 26940.00 26600.00 25600.00 26540.00 28300.00 29500.00
## 2015 27460.00 25080.00 24500.00 25200.00 24600.00 26700.00 26020.00 24960.00
## 2016 22600.00 23800.00 24700.00 25500.00 23700.00 23800.00 24300.00 26800.00
## 2017 25480.00 25000.00 26260.00 27040.00 30480.00 31780.00 30600.00 33000.00
## 2018 32520.00 30160.00 30220.00 33200.00 33120.00 35320.00 33460.00 32800.00
## 2019 33340.00 35980.00 39600.00 39400.00 37400.00 38640.00 39300.00 39580.00
## 2020 42440.00 39000.00 24100.00 25320.00 25000.00 24000.00 26800.00 26000.00
## 2021 30600.00 31400.00 28750.00 27920.00 27000.00 26400.00 27100.00 31070.00
## 2022 39500.00 39040.00 43380.00 38820.00 45320.00 33540.00 35990.00 34200.00
## 2023 41990.00 35900.00 34990.00 36000.00 27800.00 31000.00 33380.00 29100.00
## 2024 32760.00 32940.00 34280.00 33800.00 35900.00 35300.00 36520.00 37980.00
## Sep Oct Nov Dec
## 2010 29194.93 29413.54 28161.48 29314.17
## 2011 28121.73 28817.32 27426.14 28300.59
## 2012 26400.00 28660.00 27720.00 30000.00
## 2013 26900.00 26100.00 24300.00 23820.00
## 2014 27680.00 28200.00 27600.00 27640.00
## 2015 23720.00 23540.00 21660.00 20980.00
## 2016 26100.00 26500.00 24900.00 25220.00
## 2017 32740.00 28200.00 29000.00 29980.00
## 2018 31900.00 30280.00 32180.00 30400.00
## 2019 39500.00 41100.00 41680.00 44000.00
## 2020 24280.00 24500.00 27680.00 34980.00
## 2021 32960.00 33680.00 32370.00 34700.00
## 2022 31100.00 35500.00 40000.00 42500.00
## 2023 30810.00 29350.00 31720.00 33200.00
## 2024 36240.00 37800.00 38580.00 37600.00Precio de la acción mensual de Bancolombia
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("2010-01-01"), by = "month", length.out = nrow(data_col)),
y = variable1)) +
geom_line(color = "purple", 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)La gráfica muestra el precio de la acción de Bancolombia a lo largo del tiempo (2010 - 2025), expresadas en pesos colombianos.
Se puede observar del año 2015 al 2020 una tendencia alcista en el precio de la acción de Bancolombia, llegando a finales del año 2019 a un precio máximo de 44.000 (precio máximo para el periodo del 2015-2020). Posteriormente, en los primeros meses del año 2020, se presenta una caída abrupta en el precio, el cual pudo ser influenciado por el efecto de la pandemia COVID-19, llegando a un valor mínimo de 24.100. (fase de expansíón de 2 a 3 años)
Del periodo 2020 al 2025 se puede observar de igual manera una tendencia alcista, pero con muchas fluctuaciones, hasta llegar a un valor máximo histórico de 47.400, lo cual pudo ser influenciado por buenas condiciones macroeconómicas y buen desempeño del sector financiero. Cabe aclarar que este precio de la acción sobrepasa los valores de prepandemia, lo que genera buenas perspectivas hacia el año 2026.
# 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)Hay cierto componente estacional en la variable, sin embargo no es influyente, ya que los activos financieros, especialmente las acciones, no poseen estacionalidad. El residuo al ser tan grande evidencia la existencia de factores estructurales, no estacionales, como noticias que afecten al sector.
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]**SERIE ORIGINAL VS SERIE AJUSTE POR ESTACIONALIDAD
# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2010-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var1 <- ggplot() +
geom_line(aes(x = fechas_var1, y = variable1_ts), color = "purple", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Unidad de medida variable 1") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización## 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.## Warning in geom_line(aes(x = fechas_var1, y = variable1_ts), color = "purple",
## : Ignoring unknown parameters: `name`## Warning in geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", :
## Ignoring unknown parameters: `name`# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Se observa que aunque la serie ajustada sigue la serie original, en ciertos momentos hay pequeñas diferencias. Esta gráfica nos permite identificar tendencias y patrones estructurales sin la interferencia de variaciones estacionales. En conclusión, no hay mucha estacionalidad que afecte el precio de la acción y esto corrobora lo encontrado en la gráfica anterior. (hay correciones mínimas en la serie estacional).
Ahora graficamos serie original vs tendencia
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 o variable de análisis
Primero se debe obtener la tendencia de cada variable y luego graficarla
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("2010-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" = "purple", "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
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var1)La gráfica anterior permite identificar la tendencia en el precio de la acción de Bancolombia. Se observan 3 tendencias alcistas, la primera en el periodo del 2015 al 2019, luego otra del año 2015 al 2016 y la última del año 20 al 2020. De igual manera, se identificaron 2 tendencias bajistas, la primera inició en mayo del 2014, hasta abril del 2015; seguida de una tendencia bajista que inició en noviembre del 2016 hasta abril del 2018.
#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 2011
fechas_corregidas_var1 <- seq(from = as.Date("2011-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))
# Verificar longitudes
print(length(fechas_corregidas_var1))## [1] 168print(length(tasa_crecimiento_var1))## [1] 168print(length(tasa_tendencia_var1))## [1] 168Gráfico variable original y tendencia variable 1: tasa de crecimiento anual
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 = "solid") +
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)Esta gráfica nos permite observar los ciclos que sigue el precio de la acción. La tendencia en el año 2020 evidencia un ciclo contractivo, en el que hubo un decrecimiento del 30% en el precio, lo que pudo ser generado por los efectos adversos de la pandemia. La tendencia más reciente muestra que hay un ciclo expansivo, con un crecimiento del 26% anual. Tomando en cuenta que el último ciclo duró 3 años, se puede esperar que el precio de la acción vuelva a caer a finales del año 2026, por lo que se podría realizar un corto (short sell) ya que se prevee una caída en el precio.
MODEO 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 para aplicar un modelo ARIMA
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 metodología BOX-JENKINS, lo ideal es 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.
División en conjunto de entrenamiento y prueba para la variable 1 que es la elegida para pronosticar
El código siguiente divide una serie temporal (variable1_ts) en dos subconjuntos:
Conjunto de entrenamiento (train): Datos desde enero de 2012 hasta septiembre de 2024. Conjunto de prueba (test): Datos desde octubre de 2024 hasta diciembre de 2024.
Esto se hace para evaluar el desempeño de modelos de predicción en datos no vistos.
# 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(variable1_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable1_ts, end = c(2024, 9)) # Entrenamiento hasta septiembre 2024
test_ts <- window(variable1_ts, start = c(2024, 10)) # Prueba inicia desde oct2024Paso 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.
A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable 1, que es la elegida para pronosticar:
library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -3.4531, Lag order = 5, p-value = 0.04877
## alternative hypothesis: stationaryEste test arrojó un p-value de 0.04, y este valor al ser menor de 0.05 nos indica que se rechaza la hipótesisi nula, por tanto, la serie es estacionaria.
->>>>>>>Sin embargo, para repasar los temas, se hará la diferenciación.
#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 1 , una sola vez:
train_diff <- diff(train_ts, differences = 1) Diferenciación en niveles variable 1
A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable 1 para ver graficamente el cambio o ajuste:
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "purple") +
ggtitle("Variable 1:Serie Original") +
xlab("Tiempo") + ylab("Gwh")
ggplotly(p2) # Convertir en gráfico interactivo## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous. # Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], variable1_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable1_Diff)) +
geom_line(color = "red") +
ggtitle("variable 1:Serie Estacionaria (Una diferenciacion)") +
xlab("Tiempo") + ylab("variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoEjemplo 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.
A continuación se aplica la diferenciación logarítimica y la varible u objeto ahora se llama train_diff_log:
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
train_diff_log <- diff(log(train_ts), differences = 1)Ahora graficamos la serie orignal versus la serie diferenciada una vez con logaritmo
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Variable1:Serie Original") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p2) # Convertir en gráfico interactivo## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous. # 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("Variable1:Serie Estacionaria (Una diferenciacion en logaritmo)") +
xlab("Tiempo") + ylab("Variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoAhora probamos estacionariedad en la serie diferenciada ( en nivel y logaritmo)
# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_test_diff <- adf.test(train_diff)## Warning in adf.test(train_diff): p-value smaller than printed p-value print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff
## Dickey-Fuller = -6.7286, Lag order = 5, 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)## Warning in adf.test(train_diff_log): p-value smaller than printed p-value print(adf_test_diff_log)##
## Augmented Dickey-Fuller Test
##
## data: train_diff_log
## Dickey-Fuller = -6.5883, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEn el test ADF se muestra que el valor que puede tomar d=1:
El p-value ya es menor a 0.05 con una primera diferencia en ambos casos: niveles o con logaritmo natural. Por tanto el valor que puede tomar d es igual a 1
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
En el código siguiente se crean los correlogramas para determinar los posibles valores que puedeo tomar el parámetro p** y q:**
library(forecast)
# 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)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian:
p=4* q=4* (P optimo=4) (q óptimo=4)
El modelo óptimo para esta variable seria (4,1,4)
Paso 2. Estimación manual del modelo
AIC y BIC: Se usan para comparar modelos; cuanto más bajos, mejor.
Métricas de evaluación
Mean Absolute Error (MAE)
Representa el error absoluto promedio entre las predicciones del modelo y los valores reales.
Root Mean Squared Error (RMSE)
Similar al MAE, pero da más peso a los errores grandes, porque eleva las diferencias al cuadrado antes de promediarlas.
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) = 2.91%
Expresa el error en términos relativos, como porcentaje del valor real.
Interpretación: En promedio, el modelo se equivoca en un 2,91% al predecir la energía
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 2.91% sugiere un muy buen modelo para pronostico.
Estimación del modelo identificado (4,1,4)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(4,1,4)) #Se va a estimar un modelo Arima de orden (4,1,4)
summary(manual_arima_model)## Series: train_ts
## ARIMA(4,1,4)
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): Se han producido NaNs## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 0.3636 0.3758 0.7146 -0.6187 -0.4801 -0.3269 -0.8485 0.6934
## s.e. NaN NaN NaN NaN NaN NaN NaN NaN
##
## sigma^2 = 5687317: log likelihood = -1616.93
## AIC=3251.85 AICc=3252.94 BIC=3280.39
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 255.0889 2323.388 1580.52 0.4052098 5.207777 0.3382531 -0.02479758Significancia de coeficientes
library(lmtest)## Warning: package 'lmtest' was built under R version 4.4.3## Cargando paquete requerido: zoo## Warning: package 'zoo' was built under R version 4.4.3##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)## Warning in sqrt(diag(se)): Se han producido NaNs##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.36356 NaN NaN NaN
## ar2 0.37579 NaN NaN NaN
## ar3 0.71457 NaN NaN NaN
## ar4 -0.61867 NaN NaN NaN
## ma1 -0.48008 NaN NaN NaN
## ma2 -0.32689 NaN NaN NaN
## ma3 -0.84846 NaN NaN NaN
## ma4 0.69338 NaN NaN NaN##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.
- 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 demanda de energía que es la variable 1).
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 (en este caso pandemia).
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,1,4)
## Q* = 13.312, df = 16, p-value = 0.6498
##
## Model df: 8. Total lags used: 24Paso 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. En ese caso se puede mejorar el modelo, al trabajar desde el inicio con la serie ajustada por estacionalidad en caso de que se detecte un componente estacional fuerte.
Pronóstico en el test de prueba (oct, nov y dic 2024) y gráfico
#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva riable u objeto "manual_forecast"
manual_forecast <- forecast(manual_arima_model, h = length(test_ts)) #Se generan tantos pronósticos como valores tenga el conjunto de prueba (test_ts).
# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), ## Extrae las fechas del pronóstico
Pronostico = as.numeric(manual_forecast$mean), ## Valores pronosticados
Observado = as.numeric(test_ts)) ## Valores reales de la serie
# 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("Variable1:Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic2024
# 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: 2002.054cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 2048.01MAE (Mean Absolute Error)
Indica el error promedio en unidades de la variable, es decir, en número de microempresas.
RMSE (Root Mean Squared Error)
Penaliza más los errores grandes debido a la elevación al cuadrado antes de calcular la raíz.
A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic2024
# Cargar librerías necesarias
library(forecast)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## 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 37800 36115.81
## 2 2024.833 38580 36870.02
## 3 2024.917 37600 34988.00Ahora pronosticamos fuera del periodo de análisis: Enero 2025
Es decir, le sumamos al periodo de prueba (oct,nov,dic2024) una observación más (enero 2025). Es decir, se estan pronosticando 4 observaciones o meses:
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
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 36115.81
## 2 2024.833 36870.02
## 3 2024.917 34988.00
## 4 2025.000 35634.25# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_manual, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2025: 2025 = 35634.2483291662"Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual, es decir, en caso de que no se haga la dviisón inicial de conjunto de entrenamiento y prueba
Si no hay conjunto de prueba o test y solo quieres el siguiente punto u observación, el código a usar seria algo asi:
# Pronosticar octubre de 2024
future_forecast_manual <- forecast(manual_arima_model, h = 1)
# Extraer el valor específico de marzo 2025
forecast_oct2024 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para enero 2025:", forecast_oct2024))## [1] "Pronóstico para enero 2025: 36115.8140197304"“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
El modelo automático identificado es (0,1,0). Si se compara el AIC o BIC de este modelo frente el modelo manual (4,1,4), se obtiene un valor más bajo en esta métricas en este modelo automático. Probablemente pudiera ser un buen modelo para pronosticar la variable 1=energia.
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
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(0,1,0)
##
## sigma^2 = 6260946: log likelihood = -1626.92
## AIC=3255.84 AICc=3255.86 BIC=3259.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 80.23678 2495.11 1638.758 -0.03294984 5.388945 0.3507169 -0.132936Significancia de coeficientes
PONER AQUI EL CODIGO QUE ME DA ERROR
# Ajuste del modelo ARIMA(,1,0) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts,
order = c(0, 1, 0)) # Especificamos directamente (p=0, d=1, q=0)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(0,1,0)
##
## sigma^2 = 6260946: log likelihood = -1626.92
## AIC=3255.84 AICc=3255.86 BIC=3259.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 80.23678 2495.11 1638.758 -0.03294984 5.388945 0.3507169 -0.132936# 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(0,1,0)
## Q* = 24.224, df = 24, p-value = 0.4488
##
## Model df: 0. Total lags used: 24Pronóstico modelo ARIMA automático (0,1,0)
# 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("variable1")
ggplotly(p4auto) # 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_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.750 37800 36240
## 2 2024.833 38580 36240
## 3 2024.917 37600 36240POR QUÉ ME ARROJA LOS MISMOS NUMEROS EN EL PRONOSTICO????
Ahora 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_month_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_month_forecast_auto)## Tiempo Pronostico
## 1 2024.750 36240
## 2 2024.833 36240
## 3 2024.917 36240
## 4 2025.000 36240# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2025: 2025 = 36240"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el mes octubre 2024
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_auto <- future_forecast_auto$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_auto))## [1] "Pronóstico para octubre 2024: 36240"