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

install.packages("readxl")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
install.packages("tseries")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
install.packages("forecast")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
install.packages("ggplot2")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
install.packages("plotly")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
install.packages("timetk")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
library(readxl)
data_col <- read_excel("PRECIOS DEL CAFE.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric"))
head(data_col)
## # A tibble: 6 × 4
##   FECHA               PNCAFE Pinterno PECAFE
##   <dttm>               <dbl>    <dbl>  <dbl>
## 1 2012-01-01 00:00:00    535  874863.   256.
## 2 2012-02-01 00:00:00    571  826220.   246.
## 3 2012-03-01 00:00:00    576  727565.   226.
## 4 2012-04-01 00:00:00    580  703033.   215.
## 5 2012-05-01 00:00:00    689  670335.   210.
## 6 2012-06-01 00:00:00    714  592504.   186.
variable1_ts <- ts(data_col$PNCAFE, start = c(2012, 1), frequency = 12)
variable1_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012  535.0000  571.0000  576.0000  580.0000  689.0000  714.0000  668.0000
## 2013  877.0000  625.0000  617.0000  970.0000  937.0000  913.0000 1031.0000
## 2014 1011.0000  874.0000  828.0000  832.0000 1050.0000  944.0000 1236.0000
## 2015 1088.0000 1029.0000  800.0000  924.0000 1165.0000 1240.0000 1463.0000
## 2016 1136.0000 1096.0000  944.0000 1043.0000 1163.0000 1158.0000 1102.0000
## 2017 1275.0000 1293.0000 1020.0000  834.0000  901.0000 1049.0000 1373.0000
## 2018 1131.0000 1212.0000 1037.0000  874.0000 1188.0000 1087.0000 1051.0000
## 2019 1296.0000 1106.0000  914.0000 1031.0000 1115.0000 1211.0000 1317.0000
## 2020 1050.0000 1001.0000  806.0000  744.0000 1186.0000 1362.0000 1310.0000
## 2021 1081.0000 1107.0000 1050.0000  810.0000  616.0000 1052.0000 1209.0000
## 2022  868.0000  928.0000  914.0000  750.0000 1017.0000  951.0000  943.5992
## 2023  868.0000 1025.0000  799.0000  565.8672  806.1955  956.0000  947.0000
## 2024  959.0000  961.0000  865.5326  742.0000 1120.0000 1172.0000 1158.8838
##            Aug       Sep       Oct       Nov       Dec
## 2012  565.0000  519.0000  653.0000  770.0000  904.0000
## 2013  770.0000  860.0000 1058.0000 1113.0000 1115.0000
## 2014 1151.0000  912.0000 1101.0000 1115.0000 1086.0000
## 2015 1264.0000 1058.0000 1368.0000 1322.0000 1454.0000
## 2016 1189.0000 1034.0000 1395.0000 1653.0000 1319.0000
## 2017 1294.0000 1228.0000 1073.0000 1304.0000 1550.0000
## 2018 1258.0000 1050.0000 1086.0000 1300.0000 1283.0000
## 2019 1119.0000 1088.0000 1369.0000 1506.0000 1680.0000
## 2020 1091.0000  995.0000 1159.0000 1443.0000 1743.0000
## 2021  915.0000 1209.0000 1012.0000 1131.0000 1385.0000
## 2022  949.0000  834.0000  888.0000 1060.0000  981.0000
## 2023  872.0000  849.0000 1157.4585 1282.1066 1220.0000
## 2024 1049.0000 1071.2334 1339.1550 1761.4095 1798.2306
# Convertir/declarar el ISE en serie de tiempo  mensual
variable2_ts <- ts(data_col$Pinterno, start = c(2012, 1), frequency = 12)
variable2_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012  874862.9  826219.8  727564.5  703033.3  670334.7  592504.2  648096.8
## 2013  527979.8  504406.2  512520.2  515554.2  511000.0  476450.0  467209.7
## 2014  429661.3  602312.5  758746.0  796837.5  743899.2  664916.7  644649.2
## 2015  768548.4  724982.1  684923.4  687370.8  631318.5  671900.0  678778.2
## 2016  787528.2  791775.9  799129.0  756366.7  755322.6  794433.3  835516.1
## 2017  883225.8  859285.7  850064.5  804900.0  802516.1  783400.0  849322.6
## 2018  763903.2  745031.2  729854.8  715325.0  754209.7  746400.0  717838.7
## 2019  727274.2  708089.3  690580.6  680566.7  724064.5  779916.7  796483.9
## 2020  886161.3  909103.4 1143193.5 1175566.7 1068871.0  962800.0 1001451.6
## 2021 1073193.5 1113535.7 1156032.3 1207433.3 1385935.5 1420800.0 1562741.9
## 2022 2148333.3 2213333.3 1988774.2 2027448.3 2096733.3 2172233.3 2250290.3
## 2023 1838032.0 2075285.7 1993129.0 1978900.0 1895967.7 1577900.0 1318774.2
## 2024 1423000.0 1456862.0 1440129.0 1661933.3 1596129.0 1817033.3 1876387.0
##            Aug       Sep       Oct       Nov       Dec
## 2012  611621.0  623425.0  589463.7  538683.3  521262.1
## 2013  452133.1  435562.5  407205.6  384812.5  401649.2
## 2014  715621.0  715708.3  805931.5  771579.2  781746.0
## 2015  772657.3  718670.8  733637.1  735033.3  789258.1
## 2016  794032.3  860866.7  914612.9 1007533.3  860806.5
## 2017  851903.2  813762.5  776919.4  784504.2  757967.7
## 2018  705064.5  686933.3  796774.2  804283.3  727645.2
## 2019  798935.5  815450.0  819580.6  909600.0  999129.0
## 2020 1143967.7 1142233.3 1052483.9 1044700.0 1047677.4
## 2021 1704806.5 1712138.0 1784935.5 1999655.2 2116483.9
## 2022 2315548.0 2398966.7 2277290.0 1990067.0 1933032.0
## 2023 1319096.8 1285967.0 1379065.0 1406448.3 1468258.1
## 2024 1943323.0 2120233.3 2173290.0 2456566.7 2764871.0
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(data_col$PECAFE, start = c(2012, 1), frequency = 12)
variable3_ts
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2012 256.1168 246.3010 226.0661 215.2930 209.5097 185.8843 200.9219 188.5771
## 2013 168.8300 162.5343 161.7474 162.1890 160.0019 148.4153 147.5277 144.2926
## 2014 131.9742 159.6775 210.1481 216.0827 216.4171 195.8677 193.0277 210.4074
## 2015 185.3126 177.0854 155.8761 156.4597 152.3323 151.1303 145.9206 147.1319
## 2016 136.4284 137.1914 144.3697 144.5490 144.5803 153.9430 164.0368 160.7271
## 2017 163.7090 164.8057 159.2329 156.0010 151.2374 146.8607 150.9077 156.4565
## 2018 143.4810 142.3511 139.7903 138.9967 140.5035 139.1723 134.8471 130.9932
## 2019 128.6758 128.9696 125.4671 124.0727 123.1348 132.7700 138.4145 129.8942
## 2020 150.7962 144.1289 157.3050 164.6219 155.8180 149.2377 151.2900 167.6310
## 2021 175.2395 177.2779 179.1478 181.6810 196.2590 206.1064 214.4400 226.9009
## 2022 293.1955 305.5184 287.9070 292.5780 285.1871 301.2423 286.9865 292.5200
## 2023 218.2400 237.0400 228.2448 232.3695 229.0759 215.4367 190.3325 188.7070
## 2024 206.5600 210.1300 207.5470 237.3114 232.3527 249.9400 257.9800 261.0000
##           Sep      Oct      Nov      Dec
## 2012 189.2787 183.1594 171.1003 166.0471
## 2013 139.1247 135.1697 125.3610 125.9203
## 2014 208.0950 223.1000 206.7437 193.3329
## 2015 136.4220 142.7784 138.4003 139.6919
## 2016 167.6153 170.5529 179.7767 160.9187
## 2017 151.5747 145.0094 143.8483 142.5123
## 2018 126.5253 138.8045 140.7503 129.6229
## 2019 131.3550 131.8783 143.5175 160.1652
## 2020 169.8400 156.7536 161.6000 169.9900
## 2021 238.3200 257.1257 273.6352 291.9477
## 2022 296.4571 269.4900 225.1000 223.8900
## 2023 186.0000 184.9700 194.4567 206.9235
## 2024 276.2520 280.0396 293.9385 340.5167

Introduccion

Colombia es uno de los principales productores y exportadores de café en el mundo, reconocido por su café de alta calidad y producción en altura. La industria cafetera es un pilar clave de la economía colombiana, generando empleo para más de 500.000 familias y contribuyendo significativamente a las exportaciones del país. Sin embargo, el sector enfrenta desafíos como la volatilidad de los precios internacionales, los costos de producción en aumento y los efectos del cambio climático en la producción.

En este análisis, examinamos tres variables fundamentales para entender la evolución del mercado cafetero:

PNCAFE (Precio Nominal del Café - Variable 1): Representa el precio del café en dólares por quintal u otra unidad, reflejando su comportamiento en los mercados internacionales.

Pinterno (Precio Interno del Café - Variable 2): Indica el precio del café en el mercado colombiano, medido en la moneda local. Su evolución está influenciada por la oferta y demanda internas, así como por la conversión de los precios internacionales.

PECAFE (Precio de Exportación del Café - Variable 3): Muestra el valor al que se exporta el café colombiano, influenciado por acuerdos comerciales, costos logísticos y la demanda de los principales compradores internacionales.

library(ggplot2)
library(plotly)
## 
## 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("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "red", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("PN CAFE: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Serie Original

La serie original muestra la evolución del café en Colombia entre 2012 y 2024. Se observan picos en 2014 y 2021, impulsados por un incremento en la demanda global y el apoyo de la Federación Nacional de Cafeteros en la renovación de cultivos. Por otro lado, las caídas más marcadas ocurren en 2017 y 2023, influenciadas por crisis en los precios internacionales que redujeron la rentabilidad del sector. En 2017, además, factores climáticos como lluvias excesivas impactaron negativamente la producción, afectando los niveles de oferta en el mercado.

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 PN CAFE",
       x = "Tiempo",
       y = "Valor")

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

Descomposición PN CAFE

En la descomposición de la primera variable, la tendencia muestra un crecimiento constante hasta 2015, lo que refleja el impacto positivo de programas como “Renovación con Variedades Resistentes” que mejoraron la productividad. Sin embargo, a partir de 2016, la tendencia se estabiliza debido a la caída en los precios del café y el aumento en los costos de producción. La estacionalidad muestra picos recurrentes en mayo y noviembre, lo que corresponde a las cosechas principal y secundaria en Colombia. El residuo, que representa fluctuaciones no explicadas, muestra picos en 2013 y 2020, probablemente influenciados por la crisis económica de 2013 y los efectos de la pandemia de COVID-19 en 2020, que afectó la logística y exportación de café.

# 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 Pinterno del Cafe",
       x = "Tiempo",
       y = "Valor")

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

Descomposición Pinterno del Cafe

Para la segunda variable, la tendencia indica un crecimiento estable hasta 2017, seguido de una leve caída entre 2018 y 2020, lo que coincide con la reducción en la rentabilidad del sector debido a la apreciación del peso colombiano y el alza en los costos de fertilizantes. La estacionalidad presenta aumentos en abril y octubre, lo que puede deberse a la programación de exportaciones estratégicas en estos meses para aprovechar mejores precios en el mercado internacional. El residuo muestra variaciones significativas en 2008 y 2015, lo que puede estar relacionado con caídas en la producción global y restricciones climáticas en algunas regiones cafeteras del país.

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

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

Descomposición Variable 3

En esta tercera variable, la tendencia indica un crecimiento sostenido hasta 2016, año en el que Colombia alcanzó una producción récord de 14,2 millones de sacos, impulsada por inversiones en tecnología y mejora de cultivos. Sin embargo, después de ese año, la tendencia se desacelera debido a problemas estructurales como la escasez de mano de obra y la volatilidad en los precios del café. La estacionalidad muestra picos en junio y diciembre, lo que coincide con la colocación de café en mercados internacionales durante temporadas altas de consumo. El residuo muestra picos inesperados en 2012 y 2019, posiblemente relacionados con protestas de caficultores por los bajos precios y dificultades de comercialización en esos años.

# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]
# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]
# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]
# 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 = "red", 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 = "red", :
## 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.

Serie Original vs. Serie Ajustada

Al comparar la serie original con su versión ajustada, se observa que los picos más notorios en 2011, 2015 y 2021 coinciden con momentos clave para la industria cafetera en Colombia, cuando el país logró posicionarse como líder en exportaciones de café de alta calidad. En contraste, las caídas en 2009 y 2018 reflejan los efectos negativos de fenómenos climáticos y crisis de precios que redujeron la producción. La serie ajustada permite ver que, aunque hay fluctuaciones anuales, la tendencia general del sector cafetero ha sido de crecimiento, aunque con desafíos importantes en costos y rentabilidad.

# 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 = "red", 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("Variable 2:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 2") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
## Warning in geom_line(aes(x = fechas_var2, y = variable2_ts), color = "red", :
## Ignoring unknown parameters: `name`
## Warning in geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", :
## Ignoring unknown parameters: `name`
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Serie Original vs Estacionalidad Variable 2

La comparación entre la serie original y la serie ajustada del precio interno del café en Colombia muestra una tendencia de crecimiento sostenido desde el 2000 hasta 2011, con picos en 2005 y 2011, este último impulsado por problemas de oferta en Brasil y Colombia debido a condiciones climáticas adversas. A partir de 2012, la serie presenta fluctuaciones más marcadas, con descensos en 2013 y 2018, posiblemente por una mayor oferta de café en el mercado global. En 2021 y 2022, se observa un nuevo pico significativo, reflejando el impacto de la pandemia, el incremento en costos de producción y la devaluación del peso colombiano. La comparación entre la serie original y la ajustada permite ver cómo los factores económicos y climáticos han afectado el precio interno del café, mostrando una alta sensibilidad a los cambios en la oferta y demanda global.

# 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 = "red", 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("Variable 3:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 3") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
## Warning in geom_line(aes(x = fechas_var3, y = variable3_ts), color = "red", :
## Ignoring unknown parameters: `name`
## Warning in geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", :
## Ignoring unknown parameters: `name`
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

El análisis de la serie original versus la serie ajustada del precio de exportación del café en Colombia muestra variaciones significativas a lo largo del tiempo, influenciadas por factores externos como la oferta global y las políticas comerciales. Se observan picos importantes en 2011 y 2022, coincidiendo con crisis de oferta y especulación en los mercados internacionales. En 2013 y 2018, se presentan caídas notables, posiblemente debido a una mayor producción mundial, especialmente en Brasil y Vietnam, que redujo los precios internacionales. La serie ajustada suaviza las fluctuaciones de corto plazo, permitiendo identificar una tendencia general de volatilidad que afecta la competitividad del café colombiano en el mercado internacional.

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

El análisis de la serie original versus la tendencia del precio nominal del café (PNCAFE) en Colombia muestra cómo los valores observados fluctúan alrededor de una dirección general. Se destacan picos importantes en 2011 y 2022, ambos impulsados por factores externos como la oferta global restringida y la especulación en los mercados internacionales. Entre 2013 y 2018, la tendencia evidencia una caída sostenida, reflejando un exceso de oferta mundial y una menor demanda relativa. Sin embargo, desde 2020, la tendencia muestra un repunte significativo debido a problemas logísticos, costos de producción más altos y la recuperación de la demanda postpandemia. Esta comparación es clave para entender cómo los cambios estructurales afectan el mercado cafetero colombiano a largo plazo.

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("Variable 2: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 2") +
  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)

El análisis de la serie original versus la tendencia del precio interno del café (Pinterno) en Colombia refleja cómo los valores fluctúan en torno a una dirección general. Se observan incrementos significativos en 2011 y 2022, coincidiendo con picos en los precios internacionales del café y la devaluación del peso colombiano, lo que encareció el producto en el mercado interno. Entre 2013 y 2018, la tendencia muestra una caída progresiva, reflejando periodos de sobreproducción y estabilidad en la oferta local. Desde 2020, la tendencia cambia a un crecimiento sostenido debido a mayores costos de producción, inflación y un repunte en la demanda interna. Estos movimientos reflejan la relación entre la economía global y la dinámica del mercado cafetero colombiano.

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("Variable 3: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 3") +
  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)

El análisis de la serie original versus la tendencia del precio de exportación del café (PECAFE) en Colombia muestra una evolución influenciada por la oferta y demanda global. Se observan picos notables en 2011 y 2022, años en los que los precios internacionales del café alcanzaron niveles elevados debido a problemas climáticos en Brasil y otros grandes productores. Entre 2013 y 2019, la tendencia general es a la baja, reflejando una menor competitividad en los mercados internacionales y periodos de sobreoferta. A partir de 2020, la tendencia muestra un repunte significativo, impulsado por la recuperación del consumo global tras la pandemia y el encarecimiento de los costos logísticos y de producción. Estos cambios resaltan la fuerte dependencia del sector cafetero colombiano en factores externos como la especulación en mercados de futuros y la política comercial internacional.

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

El análisis de la tasa de crecimiento anual de la serie original y la tendencia del precio nominal del café (PNCAFE) en Colombia muestra fluctuaciones significativas a lo largo del tiempo. Durante los periodos de 2010-2011 y 2021-2022, la tasa de crecimiento anual alcanzó picos positivos debido a crisis de oferta en Brasil y otros países productores, lo que elevó los precios en los mercados internacionales. Sin embargo, en años como 2013-2019, la tasa de crecimiento fue negativa o cercana a cero, reflejando una caída en los precios causada por la sobreproducción y la menor demanda internacional. La tendencia general indica que, aunque hay periodos de crecimiento acelerado, la volatilidad es una característica clave del mercado cafetero, influenciada por factores como el clima, los costos de producción y la especulación en los mercados de futuros.

#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 = "red", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable2: 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 serie original y la tendencia del precio interno del café (Pinterno) en Colombia reflejan una relación estrecha con los movimientos del precio internacional y las políticas de comercialización interna. Entre 2010 y 2011, se observa un fuerte crecimiento debido al incremento en los precios internacionales, lo que impactó positivamente el valor pagado a los productores en el país. Sin embargo, entre 2013 y 2018, la tasa de crecimiento fue más moderada o negativa, evidenciando la presión de los costos de producción y la devaluación del peso colombiano, que no siempre se tradujo en mejores ingresos para los caficultores. En los años recientes (2021-2022), la tendencia muestra un repunte debido al alza global de precios, pero con una volatilidad persistente que refleja la sensibilidad del mercado interno a factores externos e internos, como subsidios, inflación y oferta local.

#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 2
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("Variable3: 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)

La tasa de crecimiento anual de la serie original y la tendencia del precio de exportación del café (PECAFE) en Colombia muestran una alta sensibilidad a los movimientos del mercado internacional. Durante el periodo 2010-2011, se observa un crecimiento significativo impulsado por el auge de los precios internacionales, alcanzando un pico histórico. Sin embargo, entre 2013 y 2019, la tasa de crecimiento muestra una tendencia decreciente debido a la sobreoferta global y la caída en la demanda en algunos mercados clave. En 2021-2022, el precio de exportación experimentó un fuerte repunte debido a la crisis en la cadena de suministros y la recuperación de la demanda global, lo que generó un impacto positivo en los ingresos de los exportadores colombianos. No obstante, la volatilidad en los años recientes sugiere que factores como la especulación en los mercados de futuros y la estabilidad de la producción nacional siguen siendo determinantes en la evolución de esta variable.

# 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)
## 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 = -3.1273, Lag order = 5, p-value = 0.1066
## 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
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")

ggplotly(acf_plot)
ggplotly(pacf_plot)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(4,2,2)) #Se va a estimar un modelo Arima de orden (1,1,1)
summary(manual_arima_model)
## Series: train_ts 
## ARIMA(4,2,2) 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2
##       0.5337  -0.2312  0.0522  -0.2214  -1.9688  0.9999
## s.e.  0.0792   0.0897  0.0892   0.0788   0.0373  0.0377
## 
## sigma^2 = 26710:  log likelihood = -986.7
## AIC=1987.4   AICc=1988.19   BIC=2008.52
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.9088543 159.1019 124.9502 -1.748565 12.08955 0.7461802
##                    ACF1
## Training set 0.02803026
library(lmtest)
## Loading required package: zoo
## 
## 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.533734   0.079217   6.7377  1.61e-11 ***
## ar2 -0.231228   0.089655  -2.5791  0.009906 ** 
## ar3  0.052237   0.089227   0.5854  0.558250    
## ar4 -0.221398   0.078822  -2.8088  0.004972 ** 
## ma1 -1.968778   0.037337 -52.7298 < 2.2e-16 ***
## ma2  0.999916   0.037748  26.4894 < 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(4,2,2)
## Q* = 83.305, df = 18, p-value = 2.249e-10
## 
## Model df: 6.   Total lags used: 24

serie temporal de residuos: Se observa una variabilidad significativa sin una tendencia clara, lo que sugiere que el modelo captura gran parte de la estructura de los datos.

Autocorrelación (ACF): Algunos lags presentan correlaciones significativas, indicando que los residuos no son completamente ruido blanco. Esto sugiere que el modelo podría mejorarse.

Histograma de residuos: La distribución se asemeja a una normal, pero con colas ligeramente pesadas. Esto indica que los residuos no son completamente gaussianos, lo que podría afectar la precisión del modelo.

#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.

El gráfico revela una discrepancia entre el pronóstico manual y el valor observado, donde el primero subestima el crecimiento real. Mientras el valor observado muestra un alza acelerada, el pronóstico crece de forma moderada. Esto sugiere la necesidad de ajustar el modelo predictivo o incorporar más factores para mejorar su precisión.

library(forecast)

# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)

# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)
## Series: train_ts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.3025  -0.4979
## s.e.   0.0786   0.0773
## 
## sigma^2 = 29051:  log likelihood = -996.22
## AIC=1998.44   AICc=1998.6   BIC=2007.51
## 
## Training set error measures:
##                    ME    RMSE      MAE        MPE    MAPE     MASE        ACF1
## Training set 14.79646 168.763 133.6028 -0.4502639 12.9936 0.797852 -0.02653252
library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.302466   0.078576 -3.8494 0.0001184 ***
## ma2 -0.497857   0.077269 -6.4432  1.17e-10 ***
## ---
## 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(4, 1, 2))  # Especificamos directamente (p=4, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(4,1,2) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ma1      ma2
##       -0.0619  -0.0201  -0.0605  -0.2810  -0.2604  -0.4371
## s.e.   0.1866   0.1323   0.0930   0.0859   0.1893   0.1699
## 
## sigma^2 = 27708:  log likelihood = -990.74
## AIC=1995.48   AICc=1996.25   BIC=2016.64
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 14.28906 162.6038 127.2701 -0.3759858 12.37715 0.7600342
##                       ACF1
## Training set -0.0005937966
# 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(4,1,2)
## Q* = 72.721, df = 18, p-value = 1.557e-08
## 
## Model df: 6.   Total lags used: 24

Serie temporal de residuos: La variabilidad de los residuos es evidente, pero sin una tendencia clara. Esto sugiere que el modelo capta bien la estructura general de los datos, aunque algunos picos extremos pueden indicar eventos no modelados.

Autocorrelación (ACF): Se observan algunos lags con correlaciones significativas, lo que indica que aún existen patrones en los residuos. Esto sugiere que el modelo ARIMA(4,1,2) no captura completamente la dinámica del proceso y podría mejorarse.

Histograma de residuos: La distribución de los residuos se asemeja a una normal, con una ligera concentración en los valores centrales. Sin embargo, las colas parecen un poco más pesadas, lo que sugiere la posible presencia de valores atípicos o fluctuaciones inesperadas.

# Generar pronóstico para el conjunto de prueba
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts))  # Predecir los valores futuros

# Crear dataframe para gráfico interactivo del pronóstico
forecast_data_auto <- data.frame(Tiempo = time(forecast_arima_auto$mean), 
                            Pronostico = as.numeric(forecast_arima_auto$mean),
                            Observado = as.numeric(test_ts))

# Graficar pronóstico junto con los valores observados reales
p4auto <- ggplot(forecast_data_auto, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("variable1")

ggplotly(p4auto)  # Convertir el gráfico en interactivo
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

El gráfico muestra una clara discrepancia entre el pronóstico y los valores observados. Mientras el pronóstico mantiene una variación leve y estable, los valores observados presentan un crecimiento acelerado. Esto indica que el modelo predictivo está subestimando la tendencia real, lo que podría deberse a una selección inadecuada de parámetros o a la omisión de factores clave que influyen en el comportamiento de la variable. Para mejorar la precisión del modelo, sería necesario recalibrarlo y considerar variables adicionales que expliquen mejor la dinámica observada.

library(forecast)
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
# 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  1339.155     1069.374
## 2 2024.833  1761.410     1031.733
## 3 2024.917  1798.231     1031.733
# 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   1069.374
## 2 2024.833   1031.733
## 3 2024.917   1031.733
## 4 2025.000   1031.733
# 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 = 1031.73275916496"
# 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: 1069.37383824435"
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(0,1,2)(0,0,2)[12] 
## 
## Coefficients:
##           ma1      ma2    sma1    sma2
##       -0.3925  -0.4243  0.3339  0.2808
## s.e.   0.0841   0.0821  0.0872  0.0688
## 
## sigma^2 = 23388:  log likelihood = -980.05
## AIC=1970.1   AICc=1970.51   BIC=1985.22
# Cargar el paquete necesario
library(forecast)

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

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(1,1,1)(1,0,2)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1    sma2
##       0.2937  -0.8450  0.9314  -0.7246  0.0519
## s.e.  0.0987   0.0495  0.0559   0.1121  0.0976
## 
## sigma^2 = 20723:  log likelihood = -972.95
## AIC=1957.9   AICc=1958.48   BIC=1976.05
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1.106225 141.1054 112.0584 -1.437845 11.04328 0.6691926
##                    ACF1
## Training set 0.03133081
# 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,1)(1,0,2)[12]
## Q* = 33.218, df = 19, p-value = 0.02268
## 
## Model df: 5.   Total lags used: 24

Gráfico de series de residuos Los residuos fluctúan alrededor de cero sin una tendencia clara lo que indica que el modelo ha logrado capturar la estructura principal de los datos Sin embargo se observan picos que reflejan momentos de alta variabilidad o posibles valores atípicos Esto sugiere que aunque el modelo es adecuado aún pueden existir patrones no completamente explicados
Función de autocorrelación (ACF) de los residuos La mayoría de los coeficientes de autocorrelación se encuentran dentro de los límites de confianza lo que indica que los residuos se comportan en su mayoría como ruido blanco No obstante algunos valores fuera de los límites sugieren que aún podría existir estructura en los datos que no ha sido completamente modelada Esto podría indicar la necesidad de ajustes adicionales en los parámetros del modelo

Histograma de los residuos con curva de densidad normal: La distribución de los residuos es aproximadamente normal con una ligera asimetría lo que indica que el supuesto de normalidad es en general válido Sin embargo la presencia de valores extremos sugiere que hay cierta dispersión que el modelo no ha capturado completamente A pesar de esto el ajuste es adecuado para la mayoría de los datos

# 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.

La gráfica muestra una clara discrepancia entre los valores pronosticados y los valores observados para la variable en cuestión. El pronóstico (línea celeste) presenta un crecimiento moderado y sostenido, mientras que los datos reales (línea naranja) evidencian un incremento más pronunciado que alcanza niveles significativamente más altos. Este contraste sugiere que el modelo de predicción podría estar subestimando la evolución de la variable, probablemente debido a la ausencia de factores o parámetros que capturen adecuadamente los cambios más bruscos en la realidad

# 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  1339.155     1186.877
## 2 2024.833  1761.410     1186.487
## 3 2024.917  1798.231     1118.730
# 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   1186.877
## 2 2024.833   1186.487
## 3 2024.917   1118.730
## 4 2025.000   1064.748

Conclusiones

Conclusiones Generales del Análisis del Sector Cafetero en Colombia (2012-2024)

Fluctuaciones en los precios del café: Se identifican ciclos de alzas y bajas en el PNCAFE (precio nominal del café), con picos en 2014 y 2021 debido a factores como el aumento de la demanda global y crisis de oferta en otros países productores. Sin embargo, las caídas en 2017 y 2023 reflejan la volatilidad del mercado internacional, afectando la rentabilidad de los caficultores.

Impacto del mercado interno: El Pinterno (precio interno del café) ha mostrado una tendencia alineada con los precios internacionales, aunque con menor volatilidad. Esto indica que las políticas de apoyo al sector, como los subsidios y compras institucionales, han amortiguado parcialmente los efectos de las caídas globales.

Comportamiento del precio de exportación: El PECAFE (precio de exportación del café) mantiene una tendencia similar a la del mercado global, reflejando la dependencia del sector colombiano en la demanda externa. Se observan periodos de crecimiento en 2014 y 2021, mientras que en 2017 y 2023 las caídas en precios limitaron los ingresos de los exportadores.

Tendencias y crecimiento: La tasa de crecimiento anual de la serie original y la tendencia muestran que el café colombiano ha mantenido un crecimiento positivo en ciertos periodos, impulsado por mejoras en productividad y estrategias de comercialización. No obstante, las fluctuaciones en la producción y los costos han generado desafíos estructurales para la industria.

Factores climáticos y económicos: Se ha identificado que eventos climáticos como el Fenómeno de La Niña (2017) han impactado negativamente la producción, reduciendo la oferta disponible en el mercado. Además, crisis económicas y devaluaciones del peso colombiano han influido en los costos de producción y en la competitividad del sector.

Recomendaciones para el sector cafetero colombiano: Diversificación de mercados para reducir la dependencia de la demanda internacional.

Inversión en tecnología y renovación de cafetales para mejorar la productividad y resiliencia climática.

Estrategias de estabilidad de precios, como coberturas financieras o acuerdos comerciales favorables.

Fortalecimiento del consumo interno, promoviendo el café colombiano en el mercado local como una alternativa rentable ante la volatilidad internacional.