library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
data_col <- read_excel("d:/Users/juans/Escritorio/Datos  Futuros ganado bovino econo.xlsx", 
    col_types = c("date", "numeric"))
# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$Último, start = c(2015, 1), frequency = 12)
variable1_ts
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2015 217.55 222.95 214.72 210.72 201.82 178.18 190.90 162.75 166.90 157.25
## 2016 145.53 147.10 144.30 140.05 140.53 123.15 121.85 128.43 130.45 122.78
## 2017 149.55 152.57 147.93 146.18 142.72 152.22 159.43 154.18 146.00 145.05
## 2018 140.18 147.40 151.32 149.32 149.45 158.18 153.47 145.22 148.85 142.10
## 2019 141.60 133.12 136.85 141.53 132.40 142.40 145.95 142.28 145.32 136.07
## 2020 119.17 135.35 132.85 144.68 140.30 141.35 137.40 141.05 138.95 137.72
## 2021 133.60 151.35 154.62 158.18 163.00 153.93 156.57 164.85 166.88 163.02
## 2022 156.35 165.13 173.60 181.55 183.48 175.68 177.63 180.48 183.70 186.15
## 2023 230.73 239.18 247.57 245.60 249.10 252.48 237.70 219.95 223.10 246.15
## 2024 244.57 258.98 259.65 256.15 240.47 244.90 245.38 259.48 262.98 275.12
##         Nov    Dec
## 2015 158.15 157.07
## 2016 125.08 133.95
## 2017 144.75 135.60
## 2018 142.88 145.25
## 2019 131.28 121.92
## 2020 138.68 143.88
## 2021 157.73 155.90
## 2022 189.80 200.82
## 2023 249.00 247.75
## 2024 273.00 286.75
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(plotly)
## Warning: package 'plotly' was built under R version 4.3.3
## 
## Attaching package: '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':
## 
##     layout
# 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("2015-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)
# 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)
# 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("2015-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
## 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 = "grey", :
## 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.
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("2015-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)
#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 marzo 2015
fechas_corregidas_var1 <- seq(from = as.Date("2015-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# Verificar longitudes
print(length(fechas_corregidas_var1))
## [1] 108
print(length(tasa_crecimiento_var1))
## [1] 108
print(length(tasa_tendencia_var1))
## [1] 108
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)
# 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
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# 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 = -2.5459, Lag order = 4, p-value = 0.3499
## alternative hypothesis: stationary
#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) 
# 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
## 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 diferenciación)") +
    xlab("Tiempo") + ylab("variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación

train_diff_log <- diff(log(train_ts), differences = 1)
# Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(Tiempo = time(train_ts), 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 diferenciación en logaritmo)") +
    xlab("Tiempo") + ylab("Variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo
# 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 = -4.7562, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
  adf_test_diff_log <- adf.test(train_diff_log)
## 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 = -4.8309, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
# 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)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,1,1)) #Se va a estimar un modelo Arima de orden (1,1,1)
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.9698  -0.9388
## s.e.  0.0545   0.0682
## 
## sigma^2 = 72.89:  log likelihood = -412.41
## AIC=830.82   AICc=831.03   BIC=839.08
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.4254525 8.427216 6.250332 0.1529181 3.767556 0.2641091
##                      ACF1
## Training set 0.0003299621
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: '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)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.969803   0.054460  17.808 < 2.2e-16 ***
## ma1 -0.938849   0.068224 -13.761 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 27.051, df = 21, p-value = 0.1692
## 
## Model df: 2.   Total lags used: 23
#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.
# 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:  12.94068
cat("RMSE Manual: ", rmse_manual, "\n")
## RMSE Manual:  13.99983
# Cargar librerías necesarias
library(forecast)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## 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
# 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    275.12     264.1888
## 2 2024.833    273.00     265.3611
## 3 2024.917    286.75     266.4980
# 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   264.1888
## 2 2024.833   265.3611
## 3 2024.917   266.4980
## 4 2025.000   267.6006
# 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 = 267.600611780108"
# 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: 264.188812919376"
# 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,2,1) 
## 
## Coefficients:
##           ma1
##       -0.9521
## s.e.   0.0291
## 
## sigma^2 = 73.05:  log likelihood = -410.6
## AIC=825.2   AICc=825.31   BIC=830.69

El modelo seleccionado por la función ARIMA es (0,2,1), lo que indica que la serie de tiempo analizada requiere dos diferenciaciones para volverse estacionaria, lo que sugiere la presencia de una tendencia fuerte en los datos originales. La ausencia de términos autoregresivos (AR = 0) implica que los valores pasados no afectan directamente las observaciones futuras, mientras que la presencia de un término de media móvil de primer orden (MA = 1) sugiere que los errores pasados tienen un impacto significativo en las predicciones actuales. El coeficiente estimado para este término MA es -0.9521, con un error estándar de 0.0291, lo que indica una fuerte correlación negativa entre los errores presentes y pasados. Los valores del Akaike Information Criterion (AIC = 825.2) y Bayesian Information Criterion (BIC = 830.69) reflejan la calidad del modelo en términos de ajuste y complejidad, mientras que el log-likelihood (-410.6) muestra el nivel de verosimilitud del modelo con los datos observados. En general, este modelo sugiere que la serie sigue una tendencia marcada y que las fluctuaciones en los datos pueden explicarse mejor a través de un modelo de media móvil simple con diferenciación de segundo orden.

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12]
darima <- Arima(train_ts, 
                order = c(0, 2, 1),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(1, 0, 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,2,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.9535  0.0475
## s.e.   0.0290  0.1093
## 
## sigma^2 = 73.58:  log likelihood = -410.51
## AIC=827.02   AICc=827.23   BIC=835.25
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1.047002 8.429769 6.304709 0.6462009 3.821528 0.2664069
##                     ACF1
## Training set -0.01899101

La parte no estacional del modelo incluye un término de media móvil (MA) de primer orden con un coeficiente de -0.9535, lo que implica que los errores pasados tienen un impacto significativo en la predicción actual. Adicionalmente, se incorpora un componente estacional con periodicidad 12, lo que sugiere la presencia de un patrón recurrente anual. En este caso, el término autorregresivo estacional (SAR(1)) tiene un coeficiente de 0.0475, indicando una leve influencia de los valores pasados dentro de la misma estación. Los valores del AIC (827.02) y BIC (835.25) muestran el balance entre ajuste y complejidad del modelo, mientras que el log-likelihood (-410.51) refleja la calidad del ajuste. En cuanto a los errores del conjunto de entrenamiento, el RMSE (8.429769) y el MAPE (3.821528%) sugieren un desempeño moderado, con cierta precisión en la predicción

# 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,2,1)(1,0,0)[12]
## Q* = 27.031, df = 21, p-value = 0.1698
## 
## Model df: 2.   Total lags used: 23

El ARIMA(0,2,1)(1,0,0)[12] sugiere que, en general, el ajuste es adecuado, aunque con algunos aspectos a considerar. En el gráfico de residuos a lo largo del tiempo, estos oscilan alrededor de cero, lo que indica que el modelo no presenta un sesgo sistemático en sus predicciones. Sin embargo, se observa cierta variabilidad en la dispersión, lo que podría sugerir la presencia de heterocedasticidad o eventos atípicos que el modelo no ha capturado completamente. En el gráfico de autocorrelación de los residuos (ACF), la mayoría de los coeficientes se encuentran dentro de los límites de significancia, lo que indica que los residuos se comportan mayormente como ruido blanco. No obstante, algunos rezagos presentan valores significativos, lo que sugiere que aún podría existir estructura en los datos no explicada por el modelo. Finalmente, el histograma de los residuos muestra una distribución aproximadamente normal, aunque con cierta asimetría y una ligera concentración de valores cercanos a cero. En conjunto, estos resultados indican que el modelo es razonablemente bueno, pero podría beneficiarse de ajustes adicionales para mejorar su capacidad predictiva.

# 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(auto_arima_model, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
  Tiempo = time(arima_forecast$mean),  # Extraer las fechas del pronóstico
  Observado = as.numeric(test_ts),  # Valores reales
  Pronosticado = as.numeric(arima_forecast$mean)  # Valores pronosticados
)

# Mostrar la tabla
print(forecast_table)
##     Tiempo Observado Pronosticado
## 1 2024.750    275.12     265.2585
## 2 2024.833    273.00     267.5369
## 3 2024.917    286.75     269.8154

Se observa que en cada punto temporal los valores observados son mayores que los pronosticados, lo que indica que el modelo de pronóstico tiende a subestimar la variable analizada. Mientras que los valores observados presentan una ligera disminución seguida de un aumento significativo, los valores pronosticados siguen una tendencia más estable y lineal. La mayor diferencia se encuentra en el último punto de tiempo, donde el valor observado alcanza 286.75, mientras que el pronosticado es de 269.8154, reflejando una subestimación de aproximadamente 16.93 unidades. Esto sugiere que el modelo de pronóstico no logra capturar completamente las fluctuaciones reales de la variable

# 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   265.2585
## 2 2024.833   267.5369
## 3 2024.917   269.8154
## 4 2025.000   272.0939
# 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 = 272.093886504647"
# 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,2,1) 
## 
## Coefficients:
##           ma1
##       -0.9521
## s.e.   0.0291
## 
## sigma^2 = 73.05:  log likelihood = -410.6
## AIC=825.2   AICc=825.31   BIC=830.69

El coeficiente del componente MA1 es de -0.9521 con un error estándar de 0.0291, lo que sugiere una fuerte dependencia negativa con los errores previos y una alta significancia estadística. La varianza del error estimada es de 73.05, mientras que la log-verosimilitud del modelo es de -410.6. Los criterios de información AIC, AICc y BIC tienen valores de 825.2, 825.31 y 830.69, respectivamente, indicando un ajuste razonable aunque no necesariamente óptimo. Dado que el modelo utiliza dos diferencias, sugiere la presencia de una fuerte tendencia en los datos originales, eliminada a través de la diferenciación.

# 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: 265.258471626162"