data <- read_excel(
  "/Users/Christian/Desktop/Universidad/Sexto Semetre/Big Data/Modulo 3/Modelo ARIMA/bogota.xlsx",
  col_types = c("date", "numeric"))
## New names:
## • `` -> `...1`
data$...1 = as.Date(data$...1)

Extración de Señales

## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

EXTRACCIÓN SEÑALES

library(ggplot2)
library(plotly)

descomp <- stl(bogts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2

dtbog <- data.frame(
  Time = rep(time(bogts), 4),  # Tiempo repetido para cada componente (son 4 componentes)
  Value = c(descomp$time.series[, "seasonal"], 
            descomp$time.series[, "trend"], 
            descomp$time.series[, "remainder"], 
            bogts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(bogts))
)


# Crear gráfico con ggplot2
p <- ggplot(dtbog, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición Variable",
       x = "Tiempo",
       y = "Valor")

ggplotly(p)

Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:

bogsa <- bogts - descomp$time.series[, "seasonal"]

Gráfico serie original VS Ajustada

# Crear vector de fechas correctamente alineado con la serie
fechabog <- seq.Date(from = as.Date("2010-01-01"), by = "month",length.out = length(bogts))

graficoaj <- ggplot() +
  geom_line(aes(x = fechabog, y = bogts), color = "darkblue", size = 0.5, linetype = "solid", name = "Serie Original" ) +
  geom_line(aes(x = fechabog, y = bogsa), color = "red", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("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))

ggplotly(graficoaj)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Análisis

Cuando se comparan las dos líneas se encuentra que la serie original y la serie ajustda por estacionalidad siguen un patrón muy similar y no existen muchas variaciones, lo que indica que el impacto estacional no es muy fuerte. Esto sugiere que las fluctuaciones en el precio de la acción están más influenciadas por eventos macroeconómicos y estructurales que por ciclos repetitivos de corto plazo.

Factores que Pueden Impactar el Comportamiento Estacional

Punto de Quiebre: Identificación y Consecuencias

En la gráfica se observan dos puntos de quiebre clave:

2020 (Pandemia de COVID-19) → Caída del 34% en cinco meses. Causa: Pánico en los mercados, menor actividad económica, reducción en créditos y aumento de la morosidad. Impacto: Bancos vieron una caída en su rentabilidad y mayor riesgo crediticio.

2021-2022 (Segunda gran caída -60%) Causa: Factores post-pandemia como políticas monetarias restrictivas, inflación alta y menor rentabilidad bancaria. Impacto: Reducción en inversiones y posible salida de capitales del mercado accionario colombiano.

¿A Quién Impacta Esta Información?

Inversionistas: Pueden ajustar su estrategia de inversión, evitando comprar en períodos de incertidumbre y aprovechando tendencias alcistas. Empresas y bancos: Podrían prepararse mejor ante ciclos de crisis y ajustar estrategias de financiamiento y préstamos.

Gobierno y reguladores: Permite tomar decisiones sobre política monetaria y medidas para evitar volatilidad excesiva en el mercado financiero.

¿Qué Se Puede Hacer con Esta Información?

Estrategias de Inversión: Identificar períodos de alta y baja volatilidad para optimizar compras y ventas de acciones.

Gestión de Riesgo: Empresas pueden anticiparse a crisis financieras y ajustar su estructura de costos y financiamiento.

Tendencia

library(ggplot2)
library(plotly)

vec <- as.numeric(bogts)
tendbbog <- as.numeric(descomp$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2010-01-01"), by = "month", length.out = length(bogts))

graficoT <- ggplot() +
  geom_line(aes(x = fechas, y = vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendbbog, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "darkblue", "Tendencia" = "black")) +
  ggtitle("Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Banco Bogota") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(graficoT)

Análisis

Proyección Futura de la Tendencia

Basándonos en la dirección de la tendencia:

Actualmente, la tendencia apunta hacia un comportamiento lateral con sesgo bajista. No se puede indentificar una recuperación en la tenencia, es decir no hay recuperación sólida en el corto plazo. Si se quiere cambiar la tendencia se necesitaria una mayor confianza de los inversionistas.

Si esta tendencia lateral se mantiene se pueden proyectar tres escenarios:

Implicaciones y Uso de Esta Información

Para inversionistas: Un mercado bajista podría ser una oportunidad para comprar si se espera una recuperación en el futuro.

Para empresas y reguladores: Importante entender qué factores pueden impulsar o frenar la recuperación del sector bancario.

Para análisis financiero: Se recomienda utilizar modelos predictivos como ARIMA o redes neuronales para evaluar escenarios futuros más precisos.

Tasa de Crecimiento

Gráfico variable original y tendencia: tasa de crecimiento anual

Análisis

La serie original muestra una evolución altamente volátil de la tasa de crecimiento anual, con grandes fluctuaciones tanto positivas como negativas.Mientras que la tendencia suaviza estos movimientos y permite identificar ciclos de crecimiento y desaceleración más estructurales.

Impacto en Diferentes Actores

Inversionistas: Una recuperación en la tasa de crecimiento puede mejorar la confianza y atraer nuevas inversiones en el Banco de Bogotá.

Empresas y Clientes: Un crecimiento estable puede indicar mejores condiciones crediticias y mayor liquidez, lo que impacta directamente a clientes y empresas que dependen de los servicios bancarios.

Reguladores y Gobierno: Una fuerte volatilidad en la tasa de crecimiento puede generar preocupaciones sobre la estabilidad del sector financiero, lo que puede llevar a regulaciones adicionales.

Proyección Futura y Uso de esta Información

Si la tendencia continúa al alza, se esperaría que la tasa de crecimiento anual se estabilice o siga mejorando en los próximos años. Sin embargo, si hay un nuevo punto de inflexión a la baja, podría ser una señal de alerta sobre posibles dificultades económicas o financieras en el sector bancario.

¿Qué se puede hacer? Para inversionistas: Evaluar el riesgo de volatilidad antes de tomar decisiones financieras. Para el banco: Implementar estrategias para mantener el crecimiento positivo, como diversificación de productos o expansión de mercados. Para reguladores: Monitorear la estabilidad del sector y prever posibles crisis financieras.

Modelo ARIMA

Dividir el conjunto de entrenamiento y prueba

Paso 1: Identificación del Modelo

Se aplica el test ADF

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -4.0233, Lag order = 5, p-value = 0.01011
## alternative hypothesis: stationary

Como el p-value es menor a 0.05 la serie es estacionaria. Por lo tanto el valor que puede tomar d es igual a 0

Identificación maunal de p y q

ggplotly(pacf_plot)

Interpretación Correlogramas

Se puede observar que los valores que pueden tomar p y q

p = 1 (Óptimo = 1)

q = 1 , 2 , 3 , 4 , 5 , 6 (Óptimo = 1)

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

Paso 2: Estimación manual del modelo

Estimación del modelo indentificado (1,0,1)

manual_arima_model <- Arima(train_ts, order = c(1,0,1))
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1       mean
##       0.9522  0.1711  60067.929
## s.e.  0.0288  0.0811   5856.177
## 
## sigma^2 = 10414845:  log likelihood = -1367.59
## AIC=2743.19   AICc=2743.47   BIC=2755.06
## 
## Training set error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 256.7795 3193.41 2327.699 0.2174524 3.779249 0.3068623 -0.03084492

**Significancia de coeficientes

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)    
## ar1       9.5224e-01 2.8814e-02 33.0472  < 2e-16 ***
## ma1       1.7112e-01 8.1071e-02  2.1108  0.03479 *  
## intercept 6.0068e+04 5.8562e+03 10.2572  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 16.955, df = 22, p-value = 0.7659
## 
## Model df: 2.   Total lags used: 24

Paso 4: Pronóstico (modelo manual)

Pronosticco en el Test de prueba

#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), 
                                   Pronostico = as.numeric(manual_forecast$mean), 
                                   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

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

A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(forecast)

# Asegurar que `manual_forecast` está generado antes de usarlo
manual_forecast <- forecast(manual_arima_model, h = length(test_ts))

# Crear `forecast_table_manual`
forecast_table_manual <- data.frame(
  Tiempo = time(manual_forecast$mean),  
  Pronostico = as.numeric(manual_forecast$mean)
)

# Convertir la columna `Tiempo` a formato Año-Mes
forecast_table_manual <- forecast_table_manual %>%
  mutate(Tiempo = as.yearmon(Tiempo))  # Convierte a "YYYY-MM"

# Mostrar la tabla corregida
print(forecast_table_manual)
##      Tiempo Pronostico
## 1  Jan 2022   69905.83
## 2  Feb 2022   69435.93
## 3  Mar 2022   68988.48
## 4  Apr 2022   68562.40
## 5  May 2022   68156.67
## 6  Jun 2022   67770.32
## 7  Jul 2022   67402.43
## 8  Aug 2022   67052.10
## 9  Sep 2022   66718.51
## 10 Oct 2022   66400.85
## 11 Nov 2022   66098.37
## 12 Dec 2022   65810.33
## 13 Jan 2023   65536.05
## 14 Feb 2023   65274.87
## 15 Mar 2023   65026.17
## 16 Apr 2023   64789.35
## 17 May 2023   64563.83
## 18 Jun 2023   64349.09
## 19 Jul 2023   64144.61
## 20 Aug 2023   63949.89
## 21 Sep 2023   63764.47
## 22 Oct 2023   63587.91
## 23 Nov 2023   63419.78
## 24 Dec 2023   63259.68
## 25 Jan 2024   63107.23
## 26 Feb 2024   62962.07
## 27 Mar 2024   62823.83
## 28 Apr 2024   62692.20
## 29 May 2024   62566.85
## 30 Jun 2024   62447.49
## 31 Jul 2024   62333.84
## 32 Aug 2024   62225.61
## 33 Sep 2024   62122.55
## 34 Oct 2024   62024.41
## 35 Nov 2024   61930.96
## 36 Dec 2024   61841.98
## 37 Jan 2025   61757.24
## 38 Feb 2025   61676.56

Ahora, ronosticamos fuera del periodo de análisis: Marzo 2025

library(forecast)
library(lubridate)
library(zoo)  # Necesario para as.yearmon

# Hacer un pronóstico para el siguiente período
next_forecast_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)

# Extraer el pronóstico y formatear la fecha como "YYYY-MM"
next_month_forecast_manual <- data.frame(
  Tiempo = format(as.Date(time(next_forecast_manual$mean) * 365.25, origin = "1970-01-01"), "%Y-%m"),  
  Pronostico = as.numeric(next_forecast_manual$mean)
)

# Mostrar el pronóstico corregido
print(next_month_forecast_manual)
##     Tiempo Pronostico
## 1  2022-01   69905.83
## 2  2022-02   69435.93
## 3  2022-03   68988.48
## 4  2022-04   68562.40
## 5  2022-05   68156.67
## 6  2022-06   67770.32
## 7  2022-07   67402.43
## 8  2022-08   67052.10
## 9  2022-09   66718.51
## 10 2022-10   66400.85
## 11 2022-11   66098.37
## 12 2022-12   65810.33
## 13 2023-01   65536.05
## 14 2023-02   65274.87
## 15 2023-03   65026.17
## 16 2023-04   64789.35
## 17 2023-05   64563.83
## 18 2023-06   64349.09
## 19 2023-07   64144.61
## 20 2023-08   63949.89
## 21 2023-09   63764.47
## 22 2023-10   63587.91
## 23 2023-11   63419.78
## 24 2023-12   63259.68
## 25 2024-01   63107.23
## 26 2024-02   62962.07
## 27 2024-03   62823.83
## 28 2024-04   62692.20
## 29 2024-05   62566.85
## 30 2024-06   62447.49
## 31 2024-07   62333.84
## 32 2024-08   62225.61
## 33 2024-09   62122.55
## 34 2024-10   62024.41
## 35 2024-11   61930.96
## 36 2024-12   61841.98
## 37 2025-01   61757.24
## 38 2025-02   61676.56
## 39 2025-03   61599.72
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_manual, 1)
print(paste("Pronóstico para Marzo de 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para Marzo de 2025: 2025-03 = 61599.7208585552"

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(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.1581
## s.e.  0.0826
## 
## sigma^2 = 10435347:  log likelihood = -1357.91
## AIC=2719.82   AICc=2719.91   BIC=2725.75
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 212.4385 3207.867 2324.834 0.3245871 3.707819 0.3064848
##                      ACF1
## Training set -0.006569839

El modelo automático identificado es (1,1,0). Si se compara el AIC o BIC de este modelo frente el modelo manual (1,0,1), se obtiene un valor más bajo en esta métricas en este modelo automático.

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1  0.15811    0.08262  1.9137  0.05566 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
darima_auto <- Arima(train_ts, 
                     order = c(1,1,0))

summary(darima_auto)
## Series: train_ts 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.1581
## s.e.  0.0826
## 
## sigma^2 = 10435347:  log likelihood = -1357.91
## AIC=2719.82   AICc=2719.91   BIC=2725.75
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 212.4385 3207.867 2324.834 0.3245871 3.707819 0.3064848
##                      ACF1
## Training set -0.006569839
checkresiduals(darima_auto)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 15.188, df = 23, p-value = 0.8878
## 
## Model df: 1.   Total lags used: 24

Pronóstico modelo ARIMA automático (1,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)
library(zoo)  # Para manejar fechas en formato "Año-Mes"

# 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 = format(as.yearmon(time(arima_forecast_auto$mean)), "%Y-%m"),  # Formato "YYYY-MM"
  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  2022-01     76800     70152.57
## 2  2022-02     72190     70145.07
## 3  2022-03     49650     70143.88
## 4  2022-04     50470     70143.69
## 5  2022-05     56920     70143.67
## 6  2022-06     46000     70143.66
## 7  2022-07     41500     70143.66
## 8  2022-08     40000     70143.66
## 9  2022-09     28540     70143.66
## 10 2022-10     26000     70143.66
## 11 2022-11     32000     70143.66
## 12 2022-12     37000     70143.66
## 13 2023-01     35420     70143.66
## 14 2023-02     25080     70143.66
## 15 2023-03     32420     70143.66
## 16 2023-04     34200     70143.66
## 17 2023-05     33800     70143.66
## 18 2023-06     32550     70143.66
## 19 2023-07     30500     70143.66
## 20 2023-08     25650     70143.66
## 21 2023-09     25200     70143.66
## 22 2023-10     24900     70143.66
## 23 2023-11     24480     70143.66
## 24 2023-12     27460     70143.66
## 25 2024-01     31800     70143.66
## 26 2024-02     32480     70143.66
## 27 2024-03     27000     70143.66
## 28 2024-04     28300     70143.66
## 29 2024-05     27540     70143.66
## 30 2024-06     26600     70143.66
## 31 2024-07     25020     70143.66
## 32 2024-08     29000     70143.66
## 33 2024-09     27720     70143.66
## 34 2024-10     28100     70143.66
## 35 2024-11     27480     70143.66
## 36 2024-12     26860     70143.66
## 37 2025-01     29500     70143.66
## 38 2025-02     30520     70143.66

Ahora pronosticamos fuera del periodo de análisis

# Cargar librerías necesarias
library(forecast)
library(lubridate)
library(zoo)
library(dplyr)

# Hacer un pronóstico para el siguiente período (trimestre, mes, etc.)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo período con formato Año-Mes
next_month_forecast_auto <- data.frame(
  Tiempo = format(as.yearmon(time(next_forecast_auto$mean)), "%Y-%m"),  # Formato YYYY-MM
  Pronostico = as.numeric(next_forecast_auto$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast_auto)
##     Tiempo Pronostico
## 1  2022-01   70152.57
## 2  2022-02   70145.07
## 3  2022-03   70143.88
## 4  2022-04   70143.69
## 5  2022-05   70143.67
## 6  2022-06   70143.66
## 7  2022-07   70143.66
## 8  2022-08   70143.66
## 9  2022-09   70143.66
## 10 2022-10   70143.66
## 11 2022-11   70143.66
## 12 2022-12   70143.66
## 13 2023-01   70143.66
## 14 2023-02   70143.66
## 15 2023-03   70143.66
## 16 2023-04   70143.66
## 17 2023-05   70143.66
## 18 2023-06   70143.66
## 19 2023-07   70143.66
## 20 2023-08   70143.66
## 21 2023-09   70143.66
## 22 2023-10   70143.66
## 23 2023-11   70143.66
## 24 2023-12   70143.66
## 25 2024-01   70143.66
## 26 2024-02   70143.66
## 27 2024-03   70143.66
## 28 2024-04   70143.66
## 29 2024-05   70143.66
## 30 2024-06   70143.66
## 31 2024-07   70143.66
## 32 2024-08   70143.66
## 33 2024-09   70143.66
## 34 2024-10   70143.66
## 35 2024-11   70143.66
## 36 2024-12   70143.66
## 37 2025-01   70143.66
## 38 2025-02   70143.66
## 39 2025-03   70143.66

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(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.1581
## s.e.  0.0826
## 
## sigma^2 = 10435347:  log likelihood = -1357.91
## AIC=2719.82   AICc=2719.91   BIC=2725.75
# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(1,0,0)
darima <- Arima(train_ts, 
                order = c(1, 1, 0),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(1, 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(1,1,0)(1,1,0)[12] 
## 
## Coefficients:
##          ar1     sar1
##       0.2147  -0.5551
## s.e.  0.0853   0.0762
## 
## sigma^2 = 13791763:  log likelihood = -1263.9
## AIC=2533.8   AICc=2533.99   BIC=2542.42
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -220.6017 3514.985 2535.068 -0.4493814 3.944172 0.3342001
##                       ACF1
## Training set -0.0004782884
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)(1,1,0)[12]
## Q* = 22.922, df = 22, p-value = 0.4061
## 
## Model df: 2.   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
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# Cargar librerías necesarias
library(forecast)
library(dplyr)
library(zoo)  # Para formatear fechas en "YYYY-MM"

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(auto_arima_model, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
  Tiempo = format(as.yearmon(time(arima_forecast$mean)), "%Y-%m"),  # Formato "YYYY-MM"
  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  2022-01     76800     70152.57
## 2  2022-02     72190     70145.07
## 3  2022-03     49650     70143.88
## 4  2022-04     50470     70143.69
## 5  2022-05     56920     70143.67
## 6  2022-06     46000     70143.66
## 7  2022-07     41500     70143.66
## 8  2022-08     40000     70143.66
## 9  2022-09     28540     70143.66
## 10 2022-10     26000     70143.66
## 11 2022-11     32000     70143.66
## 12 2022-12     37000     70143.66
## 13 2023-01     35420     70143.66
## 14 2023-02     25080     70143.66
## 15 2023-03     32420     70143.66
## 16 2023-04     34200     70143.66
## 17 2023-05     33800     70143.66
## 18 2023-06     32550     70143.66
## 19 2023-07     30500     70143.66
## 20 2023-08     25650     70143.66
## 21 2023-09     25200     70143.66
## 22 2023-10     24900     70143.66
## 23 2023-11     24480     70143.66
## 24 2023-12     27460     70143.66
## 25 2024-01     31800     70143.66
## 26 2024-02     32480     70143.66
## 27 2024-03     27000     70143.66
## 28 2024-04     28300     70143.66
## 29 2024-05     27540     70143.66
## 30 2024-06     26600     70143.66
## 31 2024-07     25020     70143.66
## 32 2024-08     29000     70143.66
## 33 2024-09     27720     70143.66
## 34 2024-10     28100     70143.66
## 35 2024-11     27480     70143.66
## 36 2024-12     26860     70143.66
## 37 2025-01     29500     70143.66
## 38 2025-02     30520     70143.66

Ahora pronosticamos fuera del periodo de análisis

# Cargar librerías necesarias
library(forecast)
library(zoo)  # Para manejar fechas en formato "Año-Mes"

# 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 con formato Año-Mes
next_month_forecast <- data.frame(
  Tiempo = format(as.yearmon(time(next_forecast$mean)), "%Y-%m"),  # Formato "YYYY-MM"
  Pronostico = as.numeric(next_forecast$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast)
##     Tiempo Pronostico
## 1  2022-01   70152.57
## 2  2022-02   70145.07
## 3  2022-03   70143.88
## 4  2022-04   70143.69
## 5  2022-05   70143.67
## 6  2022-06   70143.66
## 7  2022-07   70143.66
## 8  2022-08   70143.66
## 9  2022-09   70143.66
## 10 2022-10   70143.66
## 11 2022-11   70143.66
## 12 2022-12   70143.66
## 13 2023-01   70143.66
## 14 2023-02   70143.66
## 15 2023-03   70143.66
## 16 2023-04   70143.66
## 17 2023-05   70143.66
## 18 2023-06   70143.66
## 19 2023-07   70143.66
## 20 2023-08   70143.66
## 21 2023-09   70143.66
## 22 2023-10   70143.66
## 23 2023-11   70143.66
## 24 2023-12   70143.66
## 25 2024-01   70143.66
## 26 2024-02   70143.66
## 27 2024-03   70143.66
## 28 2024-04   70143.66
## 29 2024-05   70143.66
## 30 2024-06   70143.66
## 31 2024-07   70143.66
## 32 2024-08   70143.66
## 33 2024-09   70143.66
## 34 2024-10   70143.66
## 35 2024-11   70143.66
## 36 2024-12   70143.66
## 37 2025-01   70143.66
## 38 2025-02   70143.66
## 39 2025-03   70143.66