BANCOLOMBIA
LUISA ERAZO
2025-03-30
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.
** DESCOMPOSICIÓN DE LA SERIE TEMPORAL
# 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 principales, la primera en el periodo del 2015 al 2019, luego otra del año 2015 al 2016 y la última del año 2020 al 2022. De igual manera, se identificaron 2 tendencias bajistas, la primera inició en julio del 2019, hasta septiembre del 2020; seguida de una tendencia bajista que inició en julio del 2022 hasta septiembre del 2023.
#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] 168#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] 168#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 19.3% 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 en lo que queda del 2025 (short sell) ya que se prevee una caída en el precio.
# 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 oct2024**APLICACIÓN DE LA PRUEBA DICKEY-FULLER (ADF)
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, prima el análisis gráfico, y este demuestra que es no estacionaria, por tanto se debe diferenciar.
#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
# 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 interactivoLa gráfica permite identificar la estacionalidad de la variable, el precio de la acción se mantiene en un rango.
Ahora 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
En el código siguiente se crean los correlogramas para determinar los posibles valores que puede tomar el parámetro p** y q:**
Autocorrelation Function
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)Parcial Autocorrelation Function
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)
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.02479758Como 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) = 5.2% En este caso, un MAPE de 5.2% 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 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 el precio de la acción que es la variable 1). También se identifican picos a partir del 2023.
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 (2 rezagos), 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.
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: (aqui se tiene el pronostico de octubre del 2024, partiendo del conjunto de entrenamiento)
# 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
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.132936El 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 alto en esta métricas de este modelo automático. Lo que nos deja ver que probablemnete el modelo manual pudiera ser un buen modelo para pronosticar la variable 1=precio de la acción de Bancolombia.
# Ajuste del modelo ARIMA(0,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.El modelo automático (0,1,0) no pronostica de mejor manera dentro de prueba. No se capturan muy bien los puntos de quiebre, de hecho se observa un comportamiento lineal. No es un modelo adecuado para pornóstico fuera de muestra o a futuro.
# 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 36240Ahora 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"Modelo SARIMA automático
El modelo ajustado es un SARIMA(0,1,0). Lo que significa:
(0,1,0): Parte ARIMA no estacional: 0 términos autorregresivos (AR). 1 diferenciación (d), lo que indica que la serie fue diferenciada una vez para hacerla estacionaria. 0 términos de media móvil (MA).
# 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)
##
## sigma^2 = 6260946: log likelihood = -1626.92
## AIC=3255.84 AICc=3255.86 BIC=3259.01A continuación, se crea el objeto darima para lueg poder graficar los valores reales y observados:
# Cargar el paquete necesario
library(forecast)
# Ajustar el modelo SARIMA(0,1,0)
darima <- Arima(train_ts,
order = c(0, 1, 0), # (p,d,q) -> (0,1,1)
seasonal = list(order = c(0, 1, 0), # (P,D,Q) -> (1,0,0)
period = 12)) # Periodicidad estacional de 12 meses
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(0,1,0)(0,1,0)[12]
##
## sigma^2 = 14736669: log likelihood = -1586.19
## AIC=3174.37 AICc=3174.4 BIC=3177.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.499477 3695.175 2396.948 -0.2853013 7.784374 0.5129801
## ACF1
## Training set -0.1677151(preguntar)
# 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,1,0)(0,1,0)[12]
## Q* = 77.699, df = 24, p-value = 1.409e-07
##
## Model df: 0. Total lags used: 24Pronóstico con el modelo SARIMA
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean),
Pronostico = as.numeric(forecast_arima$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("Unidad Variable 1")
ggplotly(p4) # Convertir el gráfico en interactivo## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(forecast_arima, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
Tiempo = time(arima_forecast$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table)## Tiempo Observado Pronosticado
## 1 2024.750 37800 34780
## 2 2024.833 38580 37150
## 3 2024.917 37600 38630Ahora 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 meses.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
Tiempo = time(next_forecast$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast)## Tiempo Pronostico
## 1 2024.750 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, 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)
# 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)
##
## sigma^2 = 6260946: log likelihood = -1626.92
## AIC=3255.84 AICc=3255.86 BIC=3259.01# Pronosticar octubre 2024
future_forecast_sarima <- forecast(auto_arima_model, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_sarima <- future_forecast_sarima$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_sarima))## [1] "Pronóstico para octubre 2024: 36240"Conclusión:
El modelo manual ARIMA(4,1,4) fue el que mejor desempeño mostró en la comparación entre los datos reales y los pronosticados dentro del periodo de prueba (oct.nov.dic2024). Destacó por su mayor precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.