Introducción

En este informe, se realizará un análisis de series temporales: *1) Extracción de señales o descomposición de series temporales y 2) Pronóstico de la cantidad de toneladas de molienda de caña producidas del sector agroindustrial. El periodo de análisis es desde el primer mes de 2012 (1/2004) y el ultimo mes del 2024 (12/2024).

El análisis de series temporales permite identificar patrones subyacentes en los datos, como tendencias, estacionalidad y componentes irregulares, lo que facilita una mejor comprensión del comportamiento de la variable a analizar. Para ello, se emplearán técnicas estadísticas o métodos de descomposición y el Modelo ARIMA con el fin de extraer información clave que pueda contribuir a la toma de decisiones estratégicas.

A lo largo del informe, se presentarán los resultados obtenidos y se discutirán sus implicaciones para el ingenio Riopaila.

Instalar/Cargar librerias necesarias para el análisis

#Cargar librerías necesarias
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("C:/Extraccion de señales/Caso 2/Base Caso2 AgroIndus.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric"))
View(data_col)

Variable 1 - Exportacion de Azucares y Confiterias en USD

# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$X_AZU, start = c(2012, 1), frequency = 12)
variable1_ts
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2012 63273.73 63131.04 78703.91 59113.52 53262.61 51468.10 62038.52 80345.74
## 2013 43161.71 42983.20 42333.03 47173.48 44732.13 38836.41 47655.45 50742.94
## 2014 55511.74 58179.72 64005.94 68688.66 55205.19 55059.88 60554.66 92555.20
## 2015 49612.58 74475.97 49242.17 51675.69 43515.35 38825.29 46156.80 51083.63
## 2016 47419.21 75910.47 69847.42 35360.89 31522.37 23658.58 24070.00 44568.17
## 2017 41351.03 27715.64 42360.23 57604.69 53677.03 28307.31 49175.22 51762.84
## 2018 35806.58 58736.03 39574.41 44524.79 34029.23 37002.50 33144.20 40729.46
## 2019 40204.80 35062.81 35356.67 46959.97 34491.32 26747.84 40958.56 40604.80
## 2020 40663.19 35424.31 43779.44 34592.50 31593.03 28866.33 44149.23 42992.30
## 2021 32402.13 39991.04 54102.92 40652.57 21403.41 16218.67 39668.28 49091.77
## 2022 45770.92 42565.28 54624.16 56716.47 44781.77 48652.44 44948.21 52553.15
## 2023 54989.45 44454.87 59892.70 56043.82 64330.79 41062.23 59572.53 78671.77
## 2024 50748.42 58847.97 63213.81 65965.99 43924.39 44514.93 57855.98 51061.02
##           Sep      Oct      Nov      Dec
## 2012 87865.10 77840.99 51033.47 42670.49
## 2013 73305.51 58370.75 78884.55 81782.54
## 2014 93814.59 73510.43 89303.69 53162.63
## 2015 52296.25 56979.58 55709.22 40577.68
## 2016 44853.59 43031.82 53111.86 40480.87
## 2017 57131.74 59566.80 51835.22 48584.54
## 2018 41045.75 56862.18 57481.76 39329.39
## 2019 47657.22 61113.12 43012.01 32465.60
## 2020 53438.73 54238.00 46964.82 49451.69
## 2021 48530.61 42854.63 58422.10 56977.15
## 2022 58231.55 59529.48 61754.23 74385.71
## 2023 75900.31 84117.84 69540.10 55663.78
## 2024 72987.83 83910.01 61962.82 60231.85

Variable 2 - Produccion Molienda de Caña en Toneladas

# Convertir/declarar el ISE en serie de tiempo  mensual
variable2_ts <- ts(data_col$CAN, start = c(2012, 1), frequency = 12)
variable2_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012 1685583.9 1973654.3 2083470.2 1406867.6 1233630.5 1997317.9 2107651.9
## 2013 1415218.6 1526726.9 1484409.6 1429731.5  935326.5 1922847.6 2311142.4
## 2014 1845617.2 1984383.4 2079076.7 1833917.6 1416006.1 2065484.0 2331974.0
## 2015 1873025.6 2042447.4 2107250.0 1811593.2 1540371.6 2114224.9 2296713.9
## 2016 1883293.4 1935466.7 1949405.1 1696820.9 1325742.1 1995881.7 2264160.5
## 2017 1966713.4 2110280.0 1871801.7 1626405.6 1169223.3 2153039.2 2451640.8
## 2018 2004166.1 2083484.4 2329449.7 2001725.0 1332282.6 2212367.4 2471471.0
## 2019 1953024.0 1840400.7 2051050.6 1343556.8 1261590.0 1961549.9 2471711.0
## 2020 1983027.0 2031204.0 2157380.0 1359211.0 1238443.0 1942278.0 2198382.2
## 2021 1993003.0 2063592.0 1618398.0 1717757.0  138148.0 1926337.0 2433724.0
## 2022 2217896.7 2067395.7 1987167.2 1636843.3 1320098.2 1430592.0 2173157.6
## 2023 1673601.0 2055975.3 1770832.1 1554420.0 1001084.0 1735018.0 2138407.2
## 2024 1602184.0 1819645.0 1913302.0 1109173.0  829067.7 1721484.0 2243604.3
##            Aug       Sep       Oct       Nov       Dec
## 2012 2075408.0 2005429.3 1700439.8 1144776.2 1409399.0
## 2013 2338868.2 2292918.1 2212027.0 1871888.3 1827138.0
## 2014 2434982.1 2303703.6 2188397.9 1712510.4 2087194.9
## 2015 2339657.9 2256633.7 2184373.4 1727333.5 1917544.0
## 2016 2412379.7 2233086.3 2043861.7 1660368.1 1821465.2
## 2017 2451644.7 2327705.0 2258138.6 1834750.2 2159250.3
## 2018 2493510.6 2288295.0 2163328.9 1506036.8 2150051.0
## 2019 2484820.0 2311636.9 1997837.5 1589919.7 2065113.1
## 2020 2374975.0 2347101.9 2383201.9 1512704.7 2030826.0
## 2021 2476491.0 2354418.0 2249369.0 1758250.0 2142976.0
## 2022 2403020.5 2344196.3 2085516.8 1472515.5 1864013.8
## 2023 2142163.1 2117721.8 1673912.0 1219860.0 1785299.0
## 2024 2407235.4 2345389.6 2203425.7 1856743.9 2082758.3

Variable 3 - Produccion de Azucar en Toneladas

# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(data_col$AZUCAR, start = c(2012, 1), frequency = 12)
variable3_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012 148482.73 192196.20 202407.58 134531.82 107659.28 206133.79 222915.37
## 2013 135168.40 149319.08 142933.62 131422.29  80371.14 178655.57 230789.67
## 2014 167734.00 191996.23 200620.86 179717.02 129068.55 195426.73 237122.08
## 2015 171215.78 200101.24 202033.91 169078.05 139946.28 204938.85 237711.23
## 2016 169926.86 182905.68 179683.85 154017.36 106756.91 162952.44 203745.86
## 2017 165583.72 199198.31 173355.56 146873.54  98517.24 188802.43 222104.71
## 2018 174024.39 194809.84 230002.51 186828.20 115053.84 193215.27 232444.27
## 2019 177809.14 179759.37 198856.49 116870.98 117739.59 168612.35 226633.91
## 2020 170600.11 188936.92 202148.61 131681.81 114502.47 174447.47 203347.14
## 2021 177254.28 193573.94 145987.21 158408.21  11041.41 156211.52 220751.10
## 2022 197814.77 189215.64 180967.04 139870.52 119768.71 118137.14 189326.35
## 2023 143462.57 190188.70 162163.47 142090.67  95386.34 163133.71 198676.95
## 2024 139753.73 165380.26 184377.56 100501.17  65932.61 139433.81 185315.88
##            Aug       Sep       Oct       Nov       Dec
## 2012 226868.32 225330.03 169420.23 103433.11 138275.02
## 2013 246236.59 244437.44 232222.86 180385.01 174704.20
## 2014 263433.54 251498.08 226614.14 158069.26 198084.86
## 2015 247548.25 241929.56 224220.45 156041.36 176431.60
## 2016 229394.80 224565.44 192035.97 147874.87 156737.73
## 2017 231144.49 225560.89 213325.99 169004.01 200360.13
## 2018 243102.03 226618.82 209538.73 132131.83 197649.76
## 2019 245845.71 235865.63 202637.78 147625.70 185725.22
## 2020 228369.29 228980.14 234065.74 146351.23 193674.35
## 2021 233249.23 226637.25 211407.85 162581.72 202837.56
## 2022 217934.39 225602.05 201336.46 136489.15 178279.33
## 2023 208002.06 215285.42 163315.60 109758.97 165638.81
## 2024 219000.05 226383.47 211029.53 169151.19 189902.25

Extracción de señales

Como variables analizar en este informe se tomaron, el valor en miles dolares de las exportaciones de azucares y confiterias, la produccion en tolenadas de Molienda de caña y la produccion de azucar.

Mediante los siguientes graficos se extraeran para cada una de las variables, la estacionalidad, los errores y su tendencia, con el fin de obtener datos relevantes, que permitan a la empresa Riopaila anticiparse a los movimientos del sector cañero.

Extracción señales variable 1 ” Exportacion de Azucares y confiterias

# 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 Exportacion de Azucares",
       x = "Tiempo",
       y = "Valor")

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

La grafica de la variable 1, nos muestra las exportacion de azucares y confiteria a lo largo de los años, del 2012 al 2024, reflejando una estacionalidad, la cual alcanza su pico mas alto en los meses de Julio y su pico mas bajo en los meses de abril, para cada uno de los años analizados.Esta informacion es un dato clave para la empresa Riopaila, ya que puede permitirle anticiparse a la posible demanda de su producto.

Dado que la caña de azúcar es la principal materia prima en la producción de azúcares y confitería, a continuación se analizará el proceso de extracción de señales en la producción de caña de azúcar.

Extracción señales variable 2 Produccion de Caña de Azucar

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var2 <- stl(variable2_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var2 <- data.frame(
  Time = rep(time(variable2_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var2$time.series[, "seasonal"], 
            stl_decomp_var2$time.series[, "trend"], 
            stl_decomp_var2$time.series[, "remainder"], 
            variable2_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable2_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var2, 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 Produccion de Caña",
       x = "Tiempo",
       y = "Valor")

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

Las graficas de extracion de señales para la Variable 2 Producion de Caña, nos muestra una estacionalidad marcada en los meses de Mayo, donde alcanza su mayor nivel de produción y en marzo su nivel mas bajo.

Es importante considerar que los cultivos agrícolas siempre estarán sujetos a una carga significativa de estacionalidad, debido a factores externos como el clima, que son difíciles de controlar. Por lo tanto, los productores de este bien, como la empresa Riopapila, deben ser capaces de identificar y anticipar estos cambios a lo largo del año, con el fin de ajustar sus estrategias de manera adecuada.

Se observa en el año 2021 un pico significativo, tanto en la serie original como en la de residuos. En principio, este incremento puede atribuirse al paro nacional, que bloqueó el acceso a diversos municipios del Valle, donde se encuentran concentrados nueve de los ingenios más grandes del país.”

Extracción señales variable 3 Produccion de Azucar

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var3 <- stl(variable3_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var3 <- data.frame(
  Time = rep(time(variable3_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var3$time.series[, "seasonal"], 
            stl_decomp_var3$time.series[, "trend"], 
            stl_decomp_var3$time.series[, "remainder"], 
            variable3_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable3_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, 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 Produccion de Azucar",
       x = "Tiempo",
       y = "Valor")

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

El grafico de la variable 3 Producion de azucar, confirma la alta correlacion que existe entre la producion de caña y la producion de azucar, es decir que una buena siembra y coseha de la caña de azucar, implica una mayor produccion de sus derivados, en especial el azucar.

La grafica de estacionalidad de esta variable 3, hace evidente una estacionalidad en los meses de Mayo y marzo, meses en los que se encuentran los picos mas altos y mas bajos respectivamente de la produccion de azucar.

Si bien la tendencia suavisa algunos de los puntos mas bajos de la serie original, hace mucho mas notorios otros, como la baja producion que se dio en los primeros meses del año 2024.

Variable1 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]

Variable2 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]

Variable3 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]

Gráfico serie original VS ajustada Variable 1

# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2012-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("Exportacion de Azucar: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Miles de dólares") +
  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)

El grafico anterior, demuestra que existe una estacionalidad en la serie original que debe suavisarse, al ajustar la grafica por estacionalidad, notamos cambios significativos en algunos picos, sin embargo la grafica sigue conservando mucha simetria con la serie original. los Picos que se logran ajustar por estacionalidad se dan en el 2021,en el 2019, 2017 y el 2012.

Si bien los ingresos de exportaciones de azucares y confiteria, pueden variar segun las epocas del año y la produccion de los azucares, aun no podemos confirmar que esta variable, sea influenciada significativamente por la estacionalidad, pues la misma tambien cambia segun las demandas de los clientes y de las condiciones fluctuantes del mercado extranjero.

Gráfico serie original VS ajustada Variable 2

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
  geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Produccion de Caña: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_var2)

Esta serie ajustada, permite reducir los cambios estacionales, reducciendo significativamente la mayoria de los picos de la serie original, facilitando el análisis al eliminar la interferencia de las variaciones recurrentes, esto sin dejar a un lado los cambios estructurales importantes como el presentado en el 2021, el cual no deja de ser significativo para la precisión de las estimaciones futuras.

Gráfico serie original VS ajustada Variable 3

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var3 <- ggplot() +
  geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Produccion de Azucar: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_var3)

Como se observa en la gráfica de la producción de azúcar ajustada por estacionalidad, la extracción de señales logra ajustar significativamente la Variable 3. Aunque algunos picos persisten, como los ocurridos en 2021 debido al paro nacional, en junio de 2022 y en 2013, se debe considerar la serie ajustada para el análisis. Esto permitirá acercarse de manera más precisa a los datos reales del mercado bajo condiciones normales.

Graficas de la serie original vs tendencia

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("2012-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("Exportacion de Azucar: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Miles de dólares") +
  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)

De acuerdo con la grafica, los ingresos del sector de exportaciones de azucares y confiteria estubieron por muchos años con una tendencia a la baja, la cual se extendio desde el año 2014 hasta finalides del año 2021, donde inicio una recuperacion lenta que llego a su nivel mas alto solo hasta el año 2023. Sin embargo la tendencia aun no alcanza el pico mas alto de los ingresos historicos de estas exportaciones, el cual se dio en julio del 2014 y por lo que muestra la grafica, no se espera un crecimiento constante en los proximos meses.

si analizamos la grafica original, se observa bastante suvisada por la tendencia. esta serie cuenta con dos picos bastante bajos uno, en el año 2016 y otro en el año 2021, sin embargo, al aparecer no impactaron tan fuertemente las exportaciones de azucares.

Tendencia Variable 2

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Produccion de Caña: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_var2)

En la serie original, podemos observar una variación estacional notable, influenciada principalmente por el ciclo de producción de la caña, su recolección y las condiciones climáticas. Este comportamiento se ve interrumpido por una caída significativa en mayo de 2021, ocasionada por el paro camionero de ese año. Por otro lado, al analizar la tendencia ajustada, se logra suavizar significativamente las variaciones entre los diferentes períodos, lo que facilita la identificación de patrones reales en los datos. Al enfocarnos en esta tendencia ajustada, es posible realizar un análisis más preciso de escenarios futuros y anticipar de manera eficaz los posibles cambios o eventos que puedan impactar al sector.

Tendencia Variable 3

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Produccion de Azucar: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_var3)

Para esta gráfica, podemos observar un ajuste significativo a la estacionalidad de la producción de azúcar, lo que permite evidenciar una mayor estabilidad a lo largo del tiempo. además, se aprecia una leve tendencia a la baja entre 2022 y 2024, que podría estar influenciada por factores externos como cambios climáticos o en los costos de producción. Sin embargo, en los últimos meses se observa una recuperación, lo que sugiere un posible repunte en la producción.

Calculo de 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 2013
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# Verificar longitudes
print(length(fechas_corregidas_var1))
## [1] 144
print(length(tasa_crecimiento_var1))
## [1] 144
print(length(tasa_tendencia_var1))
## [1] 144

*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("Exportacion de Azucar: 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)

El analisis grafico de la tasa de crecimiento anual de las exportaciones de azucares y confiteria, refleja un crecimiento muy ajustado frente a la serie original. se manetiene a largo de los años una tendencia mas hacia la estabilizacion que a presentar un alza significativa.

En la serie original por el contrario vemos crecimiento bastante volatil a lo largo de los 12 años de analisis. de la imagen se puede identificar varios periodos de contracion y desaceleracion, y por lo que muestran los datos, no ha sido facil para este sector salir de ellos. la tendencia ajustada para este sector de las exportaciones muestra una leve caida al terreno negativo en el ultimo periodo de analisis, pero con una posible recuperacion, mucho mas rapida que en los periodos anteriores.

Tasa de crecimiento de la serie original vs tendencia para la variable 2

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))

# Verificar longitudes
print(length(fechas_corregidas_var2))
## [1] 144
print(length(tasa_crecimiento_var2))
## [1] 144
print(length(tasa_tendencia_var2))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Produccion de Caña: 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_var2)

La tasa de crecimiento anual de la producción de caña, en comparación con la serie original y la tendencia, muestra un crecimiento constante a lo largo de los periodos analizados, sin fluctuaciones significativas, excepto por el pico causado por el paro camionero. Aunque este evento es un cambio estructural generado por factor externo social, no refleja un impacto significativo en la tendencia de la producción de caña, pues apesar de la ocuerrencia de este evento, se observa una rápida recuperación y la tendencia general se mantiene estable.

Tasa de crecimiento de la serie original vs tendencia: variable 3

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))

# Verificar longitudes
print(length(fechas_corregidas_var3))
## [1] 144
print(length(tasa_crecimiento_var3))
## [1] 144
print(length(tasa_tendencia_var3))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 3
grafico_crecimiento_var3 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Produccion de Azucar: 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_var3)

Para la grafica de tasa de crecimiento vs la tendencia de la produccion de azucar vemos unas tasas de crecimiento muy estables tanto en la serie original como en la tendencia, se observa algunos periodos de contracion pero con impactos bajos y recuperacion instantanea, la contraccion más marcada es la del 2021, pero se debe a hechos socio enconmicos puntuales, como el paro nacional, que al final no rompe con la tendecia. el crecimiento de la produccion de azucar no nos muestra alzas significativas, a excepcion de la del 2022, pues se debe al efecto base del paro nacional del 2021, sin embargo al aislar este dato el sector se ha mantenido estable.

library(ggplot2)
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "blue", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "black", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "green", size = 0.7) +
  ggtitle("Comparativo Tasas de Crecimiento: Exportacion de Azucares, 
          Producion de Caña y Azucar") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

A partir del análisis del gráfico, se puede concluir que existe una correlación significativa entre las tres variables (exportación de azúcares, producción de caña y producción de azúcar), reflejada en la tasa de crecimiento a lo largo de los años estudiados. La relación más pronunciada se observa entre la producción de caña y la de azúcar, dado que el azúcar proviene directamente de la caña, lo que sugiere que la producción de caña influye de manera directa y proporcional en la producción de azúcar.

Asimismo, el crecimiento de las exportaciones de azúcar muestra una tendencia similar, aunque con una mayor variabilidad a lo largo de los diferentes años, lo que señala que factores externos pueden estar afectando su fluctuación, como por ejemplo el precio del azúcar a nivel internacional y otros países productores.

#Aplicacion de los modelos de prediccion

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

  1. Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro
  2. 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

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(variable2_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable2_ts, end = c(2024, 9))  # Entrenamiento hasta septiembre 2024
test_ts <- window(variable2_ts, start = c(2024, 10))  # Prueba inicia desde oct2024

Paso 1: Identificación del modelo

Identificar estacionariedad

A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable 2 Producion de Caña en Toneladas, que es la elegida para pronosticar:

library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -3.5764, Lag order = 5, p-value = 0.0378
## alternative hypothesis: stationary

El test ADF en la variable 2 arrojó un p-value igual a 0.0378, este valor es menor a 0.05, por tanto la serie si es estacionaria. sin embargo como se evidencio en la grafica de extraccion de señales, esta variable si tiene efectos de estacionalidad , por lo que se hace necesario calcular una diferenciacion 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 2

A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable 2 para ver graficamente el cambio o ajuste:

 # Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable2 = as.numeric(train_ts)), aes(x = Tiempo, y = variable2)) +
    geom_line(color = "blue") +
    ggtitle("Variable 2: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], variable2_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable2_Diff)) +
    geom_line(color = "red") +
    ggtitle("variable 2:Serie Estacionaria (Una diferenciación)") +
    xlab("Tiempo") + ylab("variable 2 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

Diferenciación en logaritmo

Aplicacion de 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)

Grafica de 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), variable2 = as.numeric(train_ts)), aes(x = Tiempo, y = variable2)) +
    geom_line(color = "blue") +
    ggtitle("Variable2:Serie Original") +
    xlab("Tiempo") + ylab("Variable2")
  
  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("Variable2:Serie Estacionaria (Una diferenciación en logaritmo)") +
    xlab("Tiempo") + ylab("Variable 2 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

Prueba de 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.8402, 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 = -7.579, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

El p-value 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

creacion de 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, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff, 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=3* 5,2,4,6 q=6*,3,2 (P optimo=3) (q óptimo=6)

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

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

Comparacion del RMSE 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) = 17.89%

Este porcentaje indica que en promedio, el modelo se equivoca en un 17,89% al predecir la produccion de Cañ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 del 17.89% sugiere un modelo aceptable para el pronostico.

Estimación del modelo identificado (3,1,6)

# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(3,1,6)) #Se va a estimar un modelo Arima de orden (3,1,6)
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(3,1,6) 
## 
## Coefficients:
##           ar1     ar2      ar3      ma1      ma2     ma3      ma4      ma5
##       -0.0005  0.0030  -0.9899  -0.5275  -0.2831  0.9149  -0.6205  -0.3403
## s.e.   0.0111  0.0125   0.0142   0.0846   0.1034  0.0719   0.0891   0.0975
##          ma6
##       0.0237
## s.e.  0.0723
## 
## sigma^2 = 8.043e+10:  log likelihood = -2124.75
## AIC=4269.5   AICc=4271.06   BIC=4299.74
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE     MASE         ACF1
## Training set 7415.553 274170.9 212547.8 -7.598575 17.8985 1.134436 -0.005749222

Significancia de coefientes

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.00052224  0.01105390  -0.0472 0.9623181    
## ar2  0.00300602  0.01246551   0.2411 0.8094410    
## ar3 -0.98988777  0.01420055 -69.7077 < 2.2e-16 ***
## ma1 -0.52745343  0.08463667  -6.2320 4.606e-10 ***
## ma2 -0.28311554  0.10337461  -2.7387 0.0061676 ** 
## ma3  0.91492049  0.07186595  12.7309 < 2.2e-16 ***
## ma4 -0.62049611  0.08911384  -6.9630 3.332e-12 ***
## ma5 -0.34033234  0.09746343  -3.4919 0.0004796 ***
## ma6  0.02369076  0.07232180   0.3276 0.7432336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación significancia coeficientes

ar3, es altamente significativo (***), lo que significa que este coeficiente tiene un impacto importante en el modelo.

Paso 3. Validación de residuos del modelo estimado manual

1. Serie de residuos (gráfico superior)

Problema posible: Se observan picos alrededor de 2021, lo que indica un cambio estructural en los datos por efecto del paro nacional.

Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)

Interpretación: La mayoría de las barras están Por fuera de las líneas azules (intervalos de confianza). lo que sugiere aún la serie tiene mucha estacionalidad.El modelo podría mejorarse pero es un modelo aceptable.

Histograma de residuos con ajuste normal (gráfico inferior derecho)

Interpretación: La curva roja representa la distribución normal teórica. La mayoria de los residuos se acercan a la normalidad, sin embargo quedan 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 el paron nacional).

Conclusión y acciones recomendadas

El modelo ARIMA (3,1,6) parece aceptable, pero tiene muchas señales que podrian mejorarse

 checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,6)
## Q* = 176.44, df = 15, p-value < 2.2e-16
## 
## Model df: 9.   Total lags used: 24

Paso 4. Pronóstico (modelo manual)

El modelo sigue la tendencia general, pero tiene un sesgo de subestimación.

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 variable 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("Variable2:Pronóstico Manual vs Observado") +
  xlab("Tiempo") + ylab("Variable2")

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

MAE (Mean Absolute Error)

Indica el error promedio en unidades de la variable, es decir, en número de toneladas de caña

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   2203426      1838684
## 2 2024.833   1856744      1340253
## 3 2024.917   2082758      1592586

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    1838684
## 2 2024.833    1340253
## 3 2024.917    1592586
## 4 2025.000    1886471
# 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 = 1886471.48400256"

Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual, es decir, en caso de que no se haga la divisó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: 1838683.5363373"

Modelo ARIMA automático

Identificación automática del modelo ARIMA

El modelo automático identificado es (1,0,2). Si se compara el AIC o BIC de este modelo frente el modelo manual (3,1,6), se obtiene un valor mas alto en esta métricas en este modelo automático. por lo anterior este no es el mejor modelo para pronosticar la variable 2 = Producion de Caña.

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,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076

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|)    
## ar1       -8.3086e-01  8.3676e-02 -9.9295 < 2.2e-16 ***
## ma1        1.4498e+00  7.8479e-02 18.4733 < 2.2e-16 ***
## ma2        6.2862e-01  7.1328e-02  8.8130 < 2.2e-16 ***
## intercept  1.9166e+06  4.5237e+04 42.3666 < 2.2e-16 ***
## ---
## 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(1, 0, 2))  # Especificamos directamente (p=4, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076
# 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(1,0,2) with non-zero mean
## Q* = 175.49, df = 21, p-value < 2.2e-16
## 
## Model df: 3.   Total lags used: 24

Pronóstico modelo ARIMA automático (1,0,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

El modelo automático (1,0,2) no pronostica mejor dentro de prueba. Hay un subestimacion sobre dos de los tres periodos analizados.

# Cargar librerías necesa
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   2203426      2096791
## 2 2024.833   1856744      1915757
## 3 2024.917   2082758      1917212

Ahora pronosticamos fuera del periodo de análisis

Es decir, le sumamos al periodo de prueba 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    2096791
## 2 2024.833    1915757
## 3 2024.917    1917212
## 4 2025.000    1916003
# 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 = 1916003.29658617"

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

Modelo SARIMA automático

Como ultima media se correo el modelo SARIMA Automatica, el cual podria ser la solución o mejora al modelo arima tradicional ya que recoge el efecto estacional de las variables, es recomendable por que la variable analizadas tienen un componente estacional fuerte.

El modelo ajustado es un SARIMA (2,0,0)(2,1,1)[12], lo que significa:

(2,0,0): Parte ARIMA no estacional: 2 términos autorregresivos (AR). sin diferenciación (d), lo que indica que la serie no fue diferenciada para hacerla estacionaria. sin término de media móvil (MA).

(2,1,1)[12]: Parte estacional con periodicidad 12 (mensual si los datos son mensuales): 2 término autorregresivo estacional (SAR). 1 diferenciaciones estacionales. 1 términos de media móvil estacionales (SMA).

# Identificación automática del mejor modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo SARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(2,0,0)(2,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1      drift
##       0.3062  0.3471  -0.2317  -0.1640  -0.6308  -198.5052
## s.e.  0.0813  0.0795   0.1466   0.1217   0.1371  1201.0995
## 
## sigma^2 = 3.591e+10:  log likelihood = -1916.58
## AIC=3847.16   AICc=3848.01   BIC=3867.81

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(2,0,0)(2,1,1)[12]
darima <- Arima(train_ts, 
                order = c(2, 0, 0),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(2, 1, 1),  # (P,D,Q) -> (1,0,0)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(2,0,0)(2,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1
##       0.3066  0.3477  -0.2323  -0.1641  -0.6300
## s.e.  0.0813  0.0794   0.1466   0.1218   0.1372
## 
## sigma^2 = 3.565e+10:  log likelihood = -1916.6
## AIC=3845.19   AICc=3845.82   BIC=3862.88
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 6199.691 178020.1 120147.6 -4.797893 11.05979 0.6412666
##                     ACF1
## Training set -0.03951852

**En el correlograma de residuos siguiente se observa que se mejora la correlación de los residuos significativamente frente a los dos modelos anteriores. Al comparar los valores reales VS pronosticados se evidencia una coinidencia mas ajustada. este modelo funciona mucho mejor que los anteriores, es por ellos que se elige como el aducuado para los pronosticos y recomendaciones para la compañia Riopaila.

# 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(2,0,0)(2,1,1)[12]
## Q* = 37.707, df = 19, p-value = 0.006466
## 
## Model df: 5.   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 2")

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(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   2203426      2139263
## 2 2024.833   1856744      1578922
## 3 2024.917   2082758      2012364

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    2139263
## 2 2024.833    1578922
## 3 2024.917    2012364
## 4 2025.000    1941393
# 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 = 1941393.1617685"

Conclusión:

El modelo automático SARIMA (2,0,0) demostró el mejor desempeño en la comparación entre los datos reales y los pronosticados durante el período de prueba (octubre, noviembre, diciembre de 2024). De los tres modelos evaluados, el SARIMA se destacó por ser el más cercano a los valores reales observados, logrando ajustarse tanto al inicio como al final del pronóstico. Este modelo fue capaz de capturar adecuadamente la tendencia general de los datos originales, incluyendo su punto de quiebre, lo que evidenció su capacidad para adaptarse a cambios estructurales en la serie temporal. Si bien el modelo SARIMA presento una ligera subestimación en los valores pronosticados, su rendimiento general superó al de los otros modelos probados, consolidándose como la opción más confiable para esta serie de datos. Estos resultados sugieren que este modelo puede ser una herramienta eficaz para futuras proyecciones, aunque sería recomendable realizar ajustes adicionales para minimizar la subestimación observada.

De acuerdo con el pronóstico calculado, se espera que la producción de caña para el mes de enero de 2025 sea de 1.941.393 toneladas, lo que representa una disminución del 6,79% en comparación con la producción de diciembre de 2024. Al analizar los datos históricos de producción, se observa que la tendencia de enero con respecto a diciembre de los últimos dos años muestra una caída en la producción, que se extienda durante el primer trimestre del año. Esta disminución es un patrón recurrente que refleja una desaceleración estacional al inicio del año. Un dato relevante para considerar es que la disminución pronosticada para enero de 2025 es aproximadamente tres puntos porcentuales inferior a la registrada en los dos últimos años, cuando la caída fue cercana al 10%. Este comportamiento podría indicar una ligera atenuación de la reducción estacional en comparación con los años previos, lo que podría ser un indicio de que factores como el clima, la gestión productiva o el rendimiento de los cultivos podrían estar contribuyendo a una menor caída en la producción al inicio del año.

Al calcular la tasa de crecimiento del pronóstico de enero de 2025 respecto a los datos reales de enero de 2024 y enero de 2023, se observa un incremento del 21,17% entre enero de 2025 y enero de 2024, y del 16% entre enero de 2025 y enero de 2023. Este comportamiento refleja una tendencia positiva para el sector agroindustrial de la caña de azúcar, ya que se anticipa que la producción de caña en enero de 2025 será superior a la de los mismos meses en los dos años anteriores. Este crecimiento no solo es alentador para el sector, sino que también tiene un impacto directo en la producción de azúcar, lo cual, a su vez, se traduce en mayores ingresos para la compañía. Este aumento en la producción podría indicar una mayor eficiencia en la cosecha, así como una recuperación o adaptación favorable frente a los desafíos del mercado o del clima, lo que refuerza la perspectiva positiva a futuro para la industria.

Ingenio Riopaila, parte del Grupo Agroindustrial Riopaila Castilla, es una de las empresas más importantes en el sector agroindustrial de Colombia, especializada en la producción de azúcar, etanol y energía a partir de la caña de azúcar. Fundada hace más de 106 años, con sede en Cali, Valle del Cauca, ha sido un actor clave en el desarrollo del sector agroindustrial colombiano y se ha consolidado como líder en la producción de azúcar, con un impacto significativo tanto a nivel nacional como internacional.

En términos de participación de mercado, Riopaila representa aproximadamente el 18% de la producción nacional de azúcar y alcohol equivalente. Se ha consolidado como uno de los principales actores del sector azucarero, su producción y comercialización se extienden a 46 países.

Riopaila no solo ha sido clave en la producción de azúcar, sino que también ha diversificado su portafolio, incursionando en la generación de energía eléctrica y la comercialización de alcohol.La empresa ha mostrado un fuerte compromiso con la sostenibilidad, implementando prácticas innovadoras como la recirculación del 97% del agua utilizada en sus procesos industriales, y protegiendo áreas naturales dentro de sus operaciones.

Los anteriores resultados del análisis de extracción de señales y de pronósticos, son relevantes para el ingenio Riopaila, en la construcción de los presupuestos de insumos, mano de obra y demás costos y gastos operativos, que impactan las operaciones del primer mes del año. Como se observó en los análisis de los datos reales históricos del mes de enero, se espera que la producción de caña en ese mes para el 2025, sea mayor que la de los años anteriores, razón por la cual, la compañía debe contar con la capacidad operativa para responder a esta mayor producción y lograr una eficiencia operativa durante este inicio del año.

Las ventas de Azúcar están directamente relacionadas con la producción de caña, actualmente el azúcar que se consume en Colombia proviene principalmente de las plantaciones nacionales que tienen los 12 ingenios del país, sin embargo, Colombia importo alrededor del 18% del azúcar que consumió internamente durante el 2024. Para Riopaila entre los muchos retos que tendrá durante este 2025, estará el de abastecer de forma eficiente y sin que le impacten significativamente los precios del azúcar del mercado internacional, que, de acuerdo con los analistas de este sector, serán bastante volátiles durante el 2025. Aunque la producción de caña disminuiría entre diciembre y enero del 2025, de acuerdo con el pronóstico, es a lo que el sector ya está acostumbrado según la data histórica, sin embargo, lo que ayudara favorablemente es la reducción de esa tasa de desaceleración en la producción, pues será la más baja de los últimos dos años. Riopaila debe prepararse para articular su estrategia operativa y la de ventas, y así buscar la manera de estabilizar los precios de su producto final, el azúcar. Es importante mencionar que podrá también verse un incremento en la producción de la energía limpia derivada del bagazo de la caña de azúcar, mercado en el que también incursiono la empresa Riopaila. la estrategia en las ventas será crucial en un contexto donde se prevé una posible mayor producción de caña en comparación con los mismos meses de años anteriores, esto puede contribuir positivamente a asegurar un suministro constante y eficiente.