Exportaciones de café Colombia
VANESSA
2025-03-19
# Cargar librerías necesarias
library(readxl) # Para leer archivos Excel
library(dplyr) # Para manipulación de datos
library(ggplot2) # Para visualización
library(plotly) # Para gráficos interactivos
library(tseries) # Para análisis de series temporales
library(forecast) # Para descomposición
library(timetk) # Para análisis y visualización de series temporalesdata_col <- read_excel("C:/Users/Nessa/OneDrive/Escritorio/UNIVERSIDAD/6 SEMESTRE/Econometría Financiera/CORTE 2/ACTIVIDADES/SECTOR_Café Mensual.xlsx")
View(data_col)# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$XCAFE, start = c(2014, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul
## 2014 139213785 199935446 149162719 219600384 208390525 157380019 225461327
## 2015 282198799 245880031 207261668 177274663 190447924 155986318 234971611
## 2016 157585765 179732103 174116820 277373927 148340613 172625931 126356056
## 2017 227673983 210303064 286061767 209098779 187682401 132054201 224405646
## 2018 247015998 182036398 195996521 197820165 184363293 141568923 158313499
## 2019 252039278 201051884 186019159 176878523 135769092 142999629 196490100
## 2020 242391717 208758804 163184893 146600070 134753520 189164759 263302402
## 2021 247990071 259822507 266853050 234504429 116711299 107484512 299415353
## 2022 320963338 370330649 399119078 315190750 262815982 352504168 379143040
## 2023 241349480 250106240 279267569 203471487 241952053 211874168 216248783
## 2024 221795367 246804767 224257517 248824091 252141311 271496247 273945661
## Aug Sep Oct Nov Dec
## 2014 213487705 212559656 247367553 261194621 239493961
## 2015 214936435 245895680 175232261 183509162 212937036
## 2016 164055478 190425586 188305907 214746415 424027561
## 2017 202232715 233587305 209847313 193015750 197817426
## 2018 205115799 172986934 186450524 186748357 209094424
## 2019 210955033 159326226 187881395 189176060 234287669
## 2020 193753299 211241673 165495693 220942574 307008770
## 2021 315159230 256436667 318369338 305963741 363128040
## 2022 308379741 323947391 312021751 281205252 336698504
## 2023 214526212 221909327 177362117 239801849 294080231
## 2024 305206809 268919784 322942527 313019803 4440550221 INTRODUCCION
El café es uno de los principales productos de exportación de Colombia, y su valor FOB (Free on Board) refleja la dinámica del comercio exterior sin costos adicionales como transporte e impuestos. Factores como la demanda global, cambios en precios internacionales y condiciones climáticas afectan su desempeño.
Sin embargo, la volatilidad de las exportaciones dificulta la identificación de tendencias claras. La extracción de señales permite aislar patrones estructurales, identificar ciclos de crecimiento y detectar eventos atípicos. Este análisis busca comprender la evolución de las exportaciones FOB de café en Colombia y sus impactos en la economía.
Con esta extraccion de señales, se espera comprender mejor los datos para hacer predicciones mas certeras en cuanto a las exportaciones en FOB de cafe en colombia.
library(ggplot2)
library(plotly)
# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable1_ts)
# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2014-01-01"), by = "month", length.out = nrow(data_col)),
y = variable1)) +
geom_line(color = "grey", linewidth = 0.4) + # Cambiado 'size' por 'linewidth'
geom_point(color = "black", size = 0.1) +
ggtitle("Variable 1: Serie original") +
xlab("Tiempo") +
ylab("Unidad Variable 1") +
theme_minimal()
ggplotly(grafico_serie)Segun esta serie original acerca de las exportaciones de cafe en Colombia, nos muestra muchas fluctuaciones a lo largo de los años. Cabe resaltar el suceso importante que todos conocemos que fue la pandemia en el año 2020. Pero segun esta grafica la caida abrupta que tuvo las exportaciones de café fue para el año 2021. Esta caida fué debido a la reduccion de la produccion nacional de café por los cambios climaticos tan adversos que enfrentaron. La produccion llegó a disminuir un 12%.
Del 2022 al 2023 se evidencia en la grafica que las exportaciones presentaba un poco problemas para su produccion y distribucion. Para el 2024-2025 ya se observa un crecimiento exponencial, en el 2024 uno de los factores fue que china se posicionó como el segundo gran comprador de café y en febrero del 2025, la produccion aumentó un 42%. Gracias a estos factores es que las exportaciones tuvieron un gran repunte.
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_var1 <- stl(variable1_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var1 <- data.frame(
Time = rep(time(variable1_ts), 4), # Tiempo repetido para cada componente (son 4 componentes)
Value = c(stl_decomp_var1$time.series[, "seasonal"],
stl_decomp_var1$time.series[, "trend"],
stl_decomp_var1$time.series[, "remainder"],
variable1_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable1_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_var1, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Descomposición temporal de la variable 1",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]2 SERIES ORIGINALES VS AJUSTADA POR ESTACIONALIDAD
# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2014-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var1 <- ggplot() +
geom_line(aes(x = fechas_var1, y = variable1_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Unidad de medida variable 1") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)3 SERIE ORIGINAL VS TENDENCIA
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("2014-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)Tasa de crecimiento de la serie original vs tendencia
#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("2014-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))
# Verificar longitudes
print(length(fechas_corregidas_var1))## [1] 120print(length(tasa_crecimiento_var1))## [1] 120print(length(tasa_tendencia_var1))## [1] 1204 TASA DE CRECIMIENTO MENSUAL: VARAIBLE ORIGINAL Y TENDENCIA VARIABLE
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)En la grafica se observa que la fuerte caida en el 2022 puede estar relacionada con opstruccion en la cadena de suministro o cambios bruscos en el mercado, aunque la tendencia nos suavisa un poco la caida. La recuperación en 2024 sugiere que el sector cafetero ha logrado adaptarse con el incremento de la produccion y volver a crecer.
Si la recuperación continúa, se podrían ver tasas de crecimiento sostenibles en los próximos años y la estabilidad en la tendencia sugiere que las exportaciones podrían estabilizarse a mediano plazo.
5 MODELO ARIMA
6 Paso 1: Identificación del modelo
Veririficamos si el modelo es estacionario o no estacionario. El test de Dickey-Fuller aumentado (ADF) se usa para verificar si una serie temporal es estacionaria, es decir, si sus propiedades estadísticas (media y varianza) permanecen constantes en el tiempo.
HO: Serie no estacionaria HI: Serie estacionaria
Si el p-valor es bajo (< 0.05) → Rechazamos la hipótesis nula y concluimos que la serie es estacionaria. Si el p-valor es alto (> 0.05) → No podemos rechazar la hipótesis nula, lo que indica que la serie NO es estacionaria.
# 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 oct2024library(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.0319, Lag order = 5, p-value = 0.1475
## alternative hypothesis: stationaryEl test ADF en la variable 1 arrojó un p-value igual a 0.14, este valor es mayor a 0.05, por tanto la serie es no estacionaria. De ese modo se debe ejecutar el código siguiente para diferenciar una vez la variable 1 y luego volver a aplicar el test ADF a esa serie diferenciada una vez:
#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 1 , una sola vez:
train_diff <- diff(train_ts, differences = 1)Diferenciación en niveles variable 1
A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable 1 para ver graficamente el cambio o ajuste:
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Variable 1: Serie Original") +
xlab("Tiempo") + ylab("Gwh")
ggplotly(p2) # Convertir en gráfico interactivo# Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], variable1_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable1_Diff)) +
geom_line(color = "red") +
ggtitle("variable 1:Serie Estacionaria (Una diferenciación)") +
xlab("Tiempo") + ylab("variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoEjemplo Diferenciación en logaritmo
la aplicación de logaritmos en series de tiempo se realiza principalmente para lograr que la serie sea más estable, lineal y estacionaria. Esta transformación es relevante porque permite modelar mejor las series que siguen un crecimiento exponencial y facilita la aplicación de técnicas estadísticas que requieren estacionariedad.
A continuación se aplica la diferenciación logarítimica y la varible u objeto ahora se llama train_diff_log:
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
train_diff_log <- diff(log(train_ts), differences = 1)# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Variable1:Serie Original") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p2) # Convertir en gráfico interactivo# Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], micro_Diff = as.numeric(train_diff_log)), aes(x = Tiempo, y = micro_Diff)) +
geom_line(color = "red") +
ggtitle("Variable1:Serie Estacionaria (Una diferenciación en logaritmo)") +
xlab("Tiempo") + ylab("Variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoAhora probamos estacionariedad en la serie diferenciada ( en nivel y logaritmo)
Ya ahora al haber corregido el modelo, se verificará si este sigue siendo estacioanrio o por el contrario sigue siendo NO estacionario
# 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 = -6.6697, 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 = -6.997, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEn el test ADF se muestra que el valor que puede tomar d=1:
El p-value ya es menor a 0.05 con una primera diferencia en ambos casos: niveles o con logaritmo natural. Por tanto el valor que puede tomar d es igual a 1
7 Identificación manual de p y q
ACF (Autocorrelation Function)
Muestra la correlación de la serie con sus rezagos.
Ayuda a determinar el parámetro q en un modelo ARIMA(p, d, q).
PACF (Partial Autocorrelation Function)
Muestra la correlación parcial entre la serie y un rezago específico, eliminando el efecto de rezagos intermedios.
Ayuda a determinar el parámetro p en un modelo ARIMA(p, d, q).
Recordemos El eje X representa los rezagos (lags). El eje Y muestra la autocorrelación parcial en cada rezago. Las líneas azules punteadas son los intervalos de confianza (aproximadamente 95%). Si una barra sobrepasa estos límites, indica una autocorrelación significativa. Si las barras caen dentro de los límites, no son significativamente diferentes de cero
En el código siguiente se crean los correlogramas para determinar los posibles valores que puedeo tomar el parámetro p y q:
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)ggplotly(pacf_plot)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian: P =
1
P = 2
Q = 1
MODELOS ARIMAS ( p, d, q): (1,1,1) (2,1,1)
8 Paso 2. Estimación manual del modelo
9 Estimación del modelo identificado (1,1,1)
library(forecast)# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,1,1)) #Se va a estimar un modelo Arima de orden (1,1,1)
summary(manual_arima_model)## Series: train_ts
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.287 -0.8106
## s.e. 0.146 0.1017
##
## sigma^2 = 2.322e+15: log likelihood = -2445.33
## AIC=4896.67 AICc=4896.86 BIC=4905.22
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3028583 47626091 36070025 -2.478215 16.72363 0.6493408
## ACF1
## Training set 0.0001801074Tenemos un MAPE del 16%, por lo que es bueno ya que no sobrepasa el 20%
Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.28696 0.14597 1.9659 0.04931 *
## ma1 -0.81064 0.10167 -7.9732 1.547e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Este modelo tiene tiene coeficientes significativos por lo que no lo vamos a descartar y continuaremos con la validacion de los residuos.
10 Paso 3. Validación de residuos del modelo estimado manual
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 34.124, df = 22, p-value = 0.04772
##
## Model df: 2. Total lags used: 24Con esta validacion de los residuos, nos muestra que el modelo esta casi en perfectas condiciones, en donde la grafica esta muy estable manteniendose sobre 0, en el AFC no muestra correlaciones significativas y el histograma se muestra mayormente centrado.
11 Paso 4. Pronóstico (modelo manual)
#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva riable u objeto "manual_forecast"
manual_forecast <- forecast(manual_arima_model, h = length(test_ts)) #Se generan tantos pronósticos como valores tenga el conjunto de prueba (test_ts).
# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), ## Extrae las fechas del pronóstico
Pronostico = as.numeric(manual_forecast$mean), ## Valores pronosticados
Observado = as.numeric(test_ts)) ## Valores reales de la serie
# Graficar pronóstico manual junto con los valores observados reales
p_manual <- ggplot(manual_forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico Manual")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Variable1:Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic 2024)
# 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: 95298151cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 112565701A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic 2024
# 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 322942527 265463685
## 2 2024.833 313019803 264471909
## 3 2024.917 444055022 264187305Ahora pronosticamos fuera del periodo de análisis: Enero 2025
# 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 265463685
## 2 2024.833 264471909
## 3 2024.917 264187305
## 4 2025.000 264105634# 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 = 264105633.627876"Aqui nos sugiere que el valor FOB de las exportaciones de cafe para el 2025 serán de $ 264.105.633 Millones de dolares
12 Modelo ARIMA automático
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)
# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)## Series: train_ts
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.5131 -0.1933
## s.e. 0.0893 0.0958
##
## sigma^2 = 2.319e+15: log likelihood = -2445.24
## AIC=4896.47 AICc=4896.66 BIC=4905.03
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 2879295 47595460 35978389 -2.526362 16.75771 0.6476911
## ACF1
## Training set -0.004734159Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.513129 0.089279 -5.7474 9.06e-09 ***
## ma2 -0.193252 0.095772 -2.0178 0.04361 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(0,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts,
order = c(0, 1, 2)) # Especificamos directamente (p=0, d=1, q=2)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.5131 -0.1933
## s.e. 0.0893 0.0958
##
## sigma^2 = 2.319e+15: log likelihood = -2445.24
## AIC=4896.47 AICc=4896.66 BIC=4905.03
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 2879295 47595460 35978389 -2.526362 16.75771 0.6476911
## ACF1
## Training set -0.004734159# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)
## Q* = 34.339, df = 22, p-value = 0.04535
##
## Model df: 2. Total lags used: 24Pronóstico modelo ARIMA automático (0,1,2)
# Generar pronóstico para el conjunto de prueba
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data_auto <- data.frame(Tiempo = time(forecast_arima_auto$mean),
Pronostico = as.numeric(forecast_arima_auto$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4auto <- ggplot(forecast_data_auto, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("variable1")
ggplotly(p4auto) # Convertir el gráfico en interactivo# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_auto <- data.frame(
Tiempo = time(arima_forecast_auto$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_auto$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_auto)## Tiempo Observado Pronosticado
## 1 2024.750 322942527 264725905
## 2 2024.833 313019803 266493567
## 3 2024.917 444055022 266493567Ahora pronosticamos fuera del periodo de análisis
# 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 264725905
## 2 2024.833 266493567
## 3 2024.917 266493567
## 4 2025.000 266493567# 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 = 266493567.375105"Para el 2025 se espera que el valor FOB de las exportaciones ocile en $266.493.567 millones de dolares
13 Modelo SARIMA automático
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts) # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)## Series: train_ts
## ARIMA(0,1,1)(0,0,2)[12]
##
## Coefficients:
## ma1 sma1 sma2
## -0.6209 0.1242 0.1711
## s.e. 0.0809 0.0898 0.0857
##
## sigma^2 = 2.273e+15: log likelihood = -2443.8
## AIC=4895.6 AICc=4895.93 BIC=4907.01A 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(0,1,1)(0,0,2)[12]
darima <- Arima(train_ts,
order = c(0, 1, 1), # (p,d,q) -> (0,1,1)
seasonal = list(order = c(0, 0, 2), # (P,D,Q) -> (0,0,2)
period = 12)) # Periodicidad estacional de 12 meses
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(0,1,1)(0,0,2)[12]
##
## Coefficients:
## ma1 sma1 sma2
## -0.6209 0.1242 0.1711
## s.e. 0.0809 0.0898 0.0857
##
## sigma^2 = 2.273e+15: log likelihood = -2443.8
## AIC=4895.6 AICc=4895.93 BIC=4907.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1900386 46928852 35469458 -2.774833 16.60656 0.6385292 0.06041693# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,0,2)[12]
## Q* = 23.429, df = 21, p-value = 0.3215
##
## Model df: 3. Total lags used: 24Pronóstico con el modelo SARIMA
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean),
Pronostico = as.numeric(forecast_arima$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("Unidad Variable 1")
ggplotly(p4) # Convertir el gráfico en interactivo# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(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 322942527 264948898
## 2 2024.833 313019803 267256623
## 3 2024.917 444055022 278910806Ahora pronosticamos fuera del periodo de análisis
# 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 264948898
## 2 2024.833 267256623
## 3 2024.917 278910806
## 4 2025.000 258183520# 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 = 258183519.682238"Para el 2025 se pronostica que el valor FOB de las exportaciones de café seran de $ 258.183.519 millones de dolares
CONCLUSION
De todos los modelos, el SARIMA (0,1,1) (0,0,2) Fue el que por lo menos intentó acercarse al real. En los modelos anteriores se observaba que la linea de pronostico se mantuvo igual en una linea recta, por el contrario en el modelo SARIMA Se nota una curva más suave y ascendente, lo que indica que el modelo está capturando mejor los patrones de la serie temporal, la diferencia entre la serie observada y el pronóstico sigue existiendo, pero parece haber una mejora en el ajuste. El modelo SARIMA mejora el pronóstico en comparación con ARIMA, porque incorpora la estacionalidad, lo cual es clave en series como las exportaciones de café que pueden tener patrones estacionales.