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) 
library(readxl)
data_col <- read_excel("ACCIONES COL.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric"))
View(data_col)
variable1_ts <- ts(data_col$BOG, start = c(2010, 1), frequency = 12)
variable1_ts
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2010 33792.66 37470.10 37072.54 39557.29 38563.39 38265.22 42042.05 49297.52
## 2011 55419.95 54167.64 54267.03 52895.42 52000.91 52935.20 51663.02 48502.39
## 2012 48661.43 50192.03 49496.29 49794.44 50032.99 49993.26 50947.37 49655.32
## 2013 55638.59 58242.64 57646.29 62615.80 69573.06 68579.19 70547.13 69553.19
## 2014 66101.69 67674.56 69645.69 69585.94 69287.31 68192.25 69685.50 72214.06
## 2015 62120.00 61880.00 55200.00 62880.00 62300.00 60000.00 60020.00 58000.00
## 2016 55800.00 56800.00 60900.00 62980.00 58900.00 58500.00 58300.00 60520.00
## 2017 59860.00 58000.00 59500.00 60300.00 62380.00 62980.00 68100.00 68980.00
## 2018 66700.00 66020.00 67920.00 69500.00 69980.00 68100.00 68500.00 69000.00
## 2019 59800.00 63120.00 66000.00 68880.00 66000.00 67220.00 71500.00 71000.00
## 2020 89000.00 86040.00 71000.00 64000.00 58000.00 60520.00 62080.00 69000.00
## 2021 83400.00 79010.00 78000.00 78500.00 70000.00 69800.00 68290.00 70500.00
## 2022 76800.00 72190.00 49650.00 50470.00 56920.00 46000.00 41500.00 40000.00
## 2023 35420.00 25080.00 32420.00 34200.00 33800.00 32550.00 30500.00 25650.00
## 2024 31800.00 32480.00 27000.00 28300.00 27540.00 26600.00 25020.00 29000.00
## 2025 29500.00 30520.00                                                      
##           Sep      Oct      Nov      Dec
## 2010 51682.89 53193.62 50987.15 57546.91
## 2011 49695.06 49357.17 47428.98 48701.16
## 2012 50410.66 54167.64 52974.94 54167.64
## 2013 67187.69 66611.25 67535.19 71178.75
## 2014 70223.06 67694.44 67000.00 66100.00
## 2015 59000.00 59000.00 58000.00 59500.00
## 2016 62600.00 62080.00 61620.00 60200.00
## 2017 68980.00 64200.00 67160.00 67460.00
## 2018 66520.00 64020.00 59400.00 55800.00
## 2019 80400.00 87020.00 86500.00 85140.00
## 2020 69000.00 68040.00 70510.00 75600.00
## 2021 72220.00 75900.00 70500.00 70200.00
## 2022 28540.00 26000.00 32000.00 37000.00
## 2023 25200.00 24900.00 24480.00 27460.00
## 2024 27720.00 28100.00 27480.00 26860.00
## 2025
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 = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Analisis Como analisis descriptivo se observa que hubo un crecimiento desde el 2010 hasta el 2022, en ese punto alcanza un pico y a inicios del 2022 empieza una caida en picada, con un intento de recuperacion en mayo, pero no logro esta recuperacion y siguio en caida, hasta ahora hay unas cuantas señales de recuperacion desde finales del 2023

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)
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]
# 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 = "grey", 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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

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

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" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:

Tasa de crecimiento de la serie de tendencia y original para la variable 1

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

# Crear vector de fechas corregido, es decir que inicie desde enero 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] 170
print(length(tasa_crecimiento_var1))
## [1] 170
print(length(tasa_tendencia_var1))
## [1] 170

Grá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 = "dashed") +
  ggtitle("Variable1: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

Analisis Se puede observar un crecimiento y una tendencia similares, la diferencia estan en los picos que alcanzo el crecimiento a diferencia de la tendencia la cual se muestra mucho mas suavizada. Observando solo los patrones de tendencia se puede ver que va a tener una pronta caida, aun asi, esa caida no va a ser significativa, en cambio esperaria un crecimiento a recuperar la constante que ha seguido la tendencia. En el 2019 el Banco de Bogota se encontraba con un crecimiento estable asi mismo como su tendencia de los años anteriores se encontraba en una constante, por eso mismo se esperaba un crecimiento constante, ademas de que logro una gran participación en BAC Credomatic, lo que consolido su presencia en Centroamerica. En el 2020 y 2021 se encontrO afectado por la pandemia (COVID-19) en el 2020 tuvo una caida debido a que este año inicio la pandemia, detuvo el consumo y afecto en general a la economia, a pesar de eso el Banco de Bogota pudo controlar la caida y no se vio tan afectado. En cambio para el 2021 se encontro con un crecimiento debido a la reactivacion de la economia, las bajas tasas de interes y el consumo. A finales del 2021 Banco de Bogota vendio su participación en BAC Credomatic. Lo cual se refleja como un decrecimiento y una tendencia a la baja, en gran parte se puede pensar que sus acciones bajaron debido a que habian perdido participacion en una entidad tan importante y ya no tenian presencia en el extranjero. En cambio se centraron en la presencia local colombiana y aumentar su liquidez. Para el 2023 se ve una recuperacion gradual, lo que indica que su estrategia ha funcionado gradualmente

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.

#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 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) - 5 # 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

#Paso 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 = -1.7511, Lag order = 5, p-value = 0.6802
## alternative hypothesis: stationary

El test ADF en la variable 1 arrojó un p-value igual a 0.68, este valor es mayor a 0.05, por tanto la serie es no estacionaria. De ese modo se debe ejecutar el código siguiente para diferenciar una vez la variable 1 y luego volver a aplicar el test ADF a esa serie diferenciada una vez:

#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 = "blue") +
    ggtitle("Variable 1:Serie Original") +
    xlab("Tiempo") + ylab("Gwh")
  
  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], variable1_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable1_Diff)) +
    geom_line(color = "red") +
    ggtitle("variable 1:Serie Estacionaria (Una diferenciación)") +
    xlab("Tiempo") + ylab("variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

##Ejemplo 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
# 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 diferenciación en logaritmo)") +
    xlab("Tiempo") + ylab("Variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

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)
  print(adf_test_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff
## Dickey-Fuller = -5.583, 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)
  print(adf_test_diff_log)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff_log
## Dickey-Fuller = -5.4138, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

En 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=2 q=2 (P optimo=1) (q óptimo=1)

El modelo óptimo para esta variable seria (2,1,2)

#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 (2,1,2)

# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(2,1,2)) #Se va a estimar un modelo Arima de orden (2,1,2)
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.5698  -0.9726  0.6118  0.9279
## s.e.   0.0405   0.0673  0.0429  0.1269
## 
## sigma^2 = 14636082:  log likelihood = -1700.15
## AIC=3410.31   AICc=3410.66   BIC=3426.16
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -34.60751 3771.291 2646.846 -0.3938437 5.110331 0.2760889
##                   ACF1
## Training set 0.1100632

Significancia de coeficientes

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 -0.569782   0.040529 -14.0585 < 2.2e-16 ***
## ar2 -0.972584   0.067268 -14.4583 < 2.2e-16 ***
## ma1  0.611841   0.042872  14.2713 < 2.2e-16 ***
## ma2  0.927866   0.126932   7.3099 2.673e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación significancia coeficientes ar1, es altamente significativos (***), lo que significa que este coeficiente tiene un impacto importante en el modelo. Como al menos uno de los dos componentes e significativo, entonces continuo con la validación de residuos del 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 observa un mayor oscile en el comportamiento desde mas o menos el 2018. esto puede ser por los cambios que han tenido desde volverse accionista mayoritario en Bac Credomatic y vender esas acciones los años siguientes.

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, hay una que sale del rango y otras estan muy cerca del borde, lo que sugiere posible estructura no capturada en los datos. 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).

Conclusión y acciones recomendadas

✅ El modelo ARIMA(1,1,1)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(2,1,2)
## Q* = 22.151, df = 20, p-value = 0.3324
## 
## Model df: 4.   Total lags used: 24

#Paso 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)

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:  2122.641
cat("RMSE Manual: ", rmse_manual, "\n")
## RMSE Manual:  2281.524

MAE (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)

# 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     28100     26272.97
## 2 2024.833     27480     28393.17
## 3 2024.917     26860     28592.48
## 4 2025.000     29500     26416.84
## 5 2025.083     30520     27462.64

Ahora 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   26272.97
## 2 2024.833   28393.17
## 3 2024.917   28592.48
## 4 2025.000   26416.84
## 5 2025.083   27462.64
## 6 2025.167   28982.75
# 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.16666666667 = 28982.7498219796"

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 octubre 2024
forecast_oct2024 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para oct 2024:", forecast_oct2024))
## [1] "Pronóstico para oct 2024: 26272.9672004904"

#Modelo ARIMA automático

##Identificación automática del modelo ARIMA

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,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.1564  -0.1524
## s.e.  0.0733   0.0678
## 
## sigma^2 = 15092385:  log likelihood = -1703.38
## AIC=3412.76   AICc=3412.9   BIC=3422.28
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -37.59856 3851.828 2653.371 -0.4281911 5.161696 0.2767696
##                      ACF1
## Training set -0.006752625

Significancia 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|)  
## ma1  0.156436   0.073324  2.1335  0.03289 *
## ma2 -0.152423   0.067794 -2.2483  0.02456 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(4,1,2) 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, 2))  # Especificamos directamente (p=0, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.1564  -0.1524
## s.e.  0.0733   0.0678
## 
## sigma^2 = 15092385:  log likelihood = -1703.38
## AIC=3412.76   AICc=3412.9   BIC=3422.28
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -37.59856 3851.828 2653.371 -0.4281911 5.161696 0.2767696
##                      ACF1
## Training set -0.006752625
# 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,2)
## Q* = 30.079, df = 22, p-value = 0.1165
## 
## Model df: 2.   Total lags used: 24

##Pronóstico modelo ARIMA automático (0,1,2)

# 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
# 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     28100     26731.00
## 2 2024.833     27480     27068.71
## 3 2024.917     26860     27068.71
## 4 2025.000     29500     27068.71
## 5 2025.083     30520     27068.71

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   26731.00
## 2 2024.833   27068.71
## 3 2024.917   27068.71
## 4 2025.000   27068.71
## 5 2025.083   27068.71
## 6 2025.167   27068.71
# 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.16666666667 = 27068.7086751677"

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: 26731.0029842982"

#Modelo SARIMA automático

# 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,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.1564  -0.1524
## s.e.  0.0733   0.0678
## 
## sigma^2 = 15092385:  log likelihood = -1703.38
## AIC=3412.76   AICc=3412.9   BIC=3422.28

A continuación, se crea el objeto darima para luego poder graficar los valores reales y observados:

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12]
darima <- Arima(train_ts, 
                order = c(0, 1, 2),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(0, 1, 2),  # (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,2)(0,1,2)[12] 
## 
## Coefficients:
##          ma1      ma2     sma1    sma2
##       0.1780  -0.1297  -0.9070  0.1002
## s.e.  0.0766   0.0716   0.0783  0.0744
## 
## sigma^2 = 16838558:  log likelihood = -1602.08
## AIC=3214.16   AICc=3214.54   BIC=3229.66
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -442.4427 3901.448 2695.434 -1.067584 5.147938 0.2811572
##                     ACF1
## Training set -0.01745897
# 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,2)(0,1,2)[12]
## Q* = 30.285, df = 20, p-value = 0.06536
## 
## Model df: 4.   Total lags used: 24

Pronó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
# 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     28100     27149.80
## 2 2024.833     27480     27431.62
## 3 2024.917     26860     29042.78
## 4 2025.000     29500     31419.04
## 5 2025.083     30520     28299.62

Ahora 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   26731.00
## 2 2024.833   27068.71
## 3 2024.917   27068.71
## 4 2025.000   27068.71
## 5 2025.083   27068.71
## 6 2025.167   27068.71
# 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.16666666667 = 27068.7086751677"

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,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.1564  -0.1524
## s.e.  0.0733   0.0678
## 
## sigma^2 = 15092385:  log likelihood = -1703.38
## AIC=3412.76   AICc=3412.9   BIC=3422.28
# 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: 26731.0029842982"

Conclusión: Ningun modelo ha mostrado una semejanza al pronostico observado, por ende se concluye que ningun modelo funciona y es mejor realizar un modelo multivariado