Instalar/Cargar librerias necesarias para el análisis

#Cargar librerías necesarias
library(readxl)  
library(tseries) 
library(forecast) 
library(ggplot2)  
library(plotly)  
library(timetk)

Introducción

Para efectos del siguiente análisis económico y financiero a partir de la extracción de señales y pronóstico ARIMA, me centraré en la revisión del impacto de cuatro indicadores claves en el sector productivo del cartón y el papel, dado que, es una de las industrias claves del valle del cauca y del país, ya que contribuye no solamente a la innovación, sostenibilidad ambienta y generación de empleo. Este sector contribuye con un porcentaje del 4,6% en PIB nacional y si bien, sus exportaciones están enfocadas generalmente en sur américa y Centroamérica, mediante este análisis podremos identificar varios factores como: cómo está la producción nacional, qué de esa producción nacional se consume en el interior del país y cual se exporta, correlación del consumo de energía en el sector manufacturero en el valle frente a la producción de cartón y a su vez, la producción en específico del cartón en el valle. De conformidad con lo informado por diarios como el País, el Valle del Cauca concentra entre el 70 % y 80 % de la producción de papel de Colombia, siendo el principal departamento productor de papel y a su vez, también se ve impactado porque por el puerto de Buenaventura se realizan las exportaciones de este insumo generando un impacto directo en las exportaciones, con cifras como US$163,2 millones en ventas de papel y cartón, representando alrededor del 4,1 % de las exportaciones regionales.

Metodología

Cargar base de datos

library(readxl)
Base_Caso2_Trabajo <- read_excel("C:/Users/Jéssica Torres L/OneDrive/Documentos/Maestría en Finanzas/Entrega 2/Base Caso2 Trabajo.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric"))

PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):

Variable 1: ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

# Convertir/declarar variable 1=ENER_v en serie de tiempo mensual
variable1_ts <- ts(Base_Caso2_Trabajo$ENER_V, start = c(2012, 1), frequency = 12)

Variable 2:COR_V Empaques de cartón en el Valle

# Convertir/declarar el COR_V en serie de tiempo  mensual
variable2_ts <- ts(Base_Caso2_Trabajo$COR_V, start = c(2012, 1), frequency = 12)

Variable 3: X_PAPEL Indicador de Exportaciones de papel y cartón

# Convertir/declarar las exportaciones de papel y cartón en serie de tiempo mensual
variable3_ts <- ts(Base_Caso2_Trabajo$X_PAPEL, start = c(2012, 1), frequency = 12)

Variable 4: IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia

# Convertir/declarar Indice de producción industrial regional de papel y cartón en serie de tiempo mensual
variable4_ts <- ts(Base_Caso2_Trabajo$IPIR_PAPEL, start = c(2018, 1), frequency = 12)

Extracción de señales

Gráfico inicial de la variable 1: ENER_V (Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle en niveles) -Original

library(ggplot2)
library(plotly)

# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
Base_Caso2_Trabajo$variable1 <- as.numeric(variable1_ts)

# Crear el gráfico
grafico_serie <- ggplot(Base_Caso2_Trabajo, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(Base_Caso2_Trabajo)), 
                                      y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: ENER_V Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales variable 1: ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

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

# Manejo de NA
library(zoo)
variable1_ts <- na.approx(variable1_ts)

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

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

Extracción señales variable 2: COR_V Empaques de cartón en el Valle

# 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 variable 2: COR_v",
       x = "Tiempo",
       y = "Valor")

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

Extracción señales variable 3: X_PAPEL Indicador de Exportaciones de papel y cartón

# 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 variable 3: X_PAPEL",
       x = "Tiempo",
       y = "Valor")

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

Extracción señales variable 4: IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia

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

# Manejo de NA
library(zoo)
variable4_ts <- na.approx(variable1_ts)

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

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

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var4, 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 4: IPIR_PAPEL ",
       x = "Tiempo",
       y = "Valor")

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

Variables ajustadas por estacionalidad

Variable 1: ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle, ajustada por estacionalidad

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

Variable 2:COR_V Empaques de cartón en el Valle, ajustada por estacionalidad

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

Variable 3: X_PAPEL Indicador de Exportaciones de papel y cartón, ajustada por estacionalidad

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

Variable 4: IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia, ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable4_sa <- variable4_ts - stl_decomp_var4$time.series[, "seasonal"]

Ahora si se puede graficar las series originales versus la ajustada por estacionalidad

Gráfico serie original VS ajustada Variable 1:ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

# 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("Variable 1 ENER_V: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)

Gráfico serie original VS ajustada Variable 2:COR_V Empaques de cartón en el Valle

# 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("Variable 2 COR_V: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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

Gráfico serie original VS ajustada Variable 3:X_PAPEL Indicador de Exportaciones de papel y cartón

# 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("Variable 3 X_PAPEL: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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

Gráfico serie original VS ajustada Variable 4:IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var4 <- ggplot() +
  geom_line(aes(x = fechas_var4, y = variable4_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var4, y = variable4_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 4 IPIR_PAPEL:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 4") +
  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_var4)

Ahora, serie original vs tendencia

Primero se debe obtener la tendencia de cada variable y luego graficarla

Tendencia Variable 1: ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

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 ENER_V: 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)

Tendencia Variable 2:COR_V Empaques de cartón en el Valle

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 COR_V: 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)

Tendencia Variable 3: X_PAPEL Indicador de Exportaciones de papel y cartón

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 X_PAPEL: 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)

Tendencia Variable 4: IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable4_vec <- as.numeric(variable4_ts)
tendencia_var4 <- as.numeric(stl_decomp_var4$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var4 <- ggplot() +
  geom_line(aes(x = fechas, y = variable4_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var4, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 4 IPIR_PAPEL: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 4") +
  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_var4)

Tasa de crecimiento de la serie de tendencia y original para la variable 1, ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

#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: ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle - tasa de crecimiento anual**

library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1 ENER_V
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 ENER_V: 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)

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2: COR_V Empaques de cartón en el Valle

#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: COR_v
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("Variable2 COR_v: 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)

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3 X_PAPEL Indicador de Exportaciones de papel y cartón

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

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 4 IPIR_PAPEL Índice de Producción Industrial Regional (papel y cartón) de Colombia

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var4 <- (variable4_ts[(13:length(variable4_ts))] / variable4_ts[1:(length(variable4_ts) - 12)] - 1) * 100
tasa_tendencia_var4 <- (tendencia_var4[(13:length(tendencia_var4))] / tendencia_var4[1:(length(tendencia_var4) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var4 <- seq(from = as.Date("2019-01-01"), by = "month", length.out = length(tasa_crecimiento_var4))

# Verificar longitudes
print(length(fechas_corregidas_var4))

## [1] 144

print(length(tasa_crecimiento_var4))

## [1] 144

print(length(tasa_tendencia_var4))

## [1] 144

# Gráfico de la tasa de crecimiento anual variable 4 IPIR_PAPEL
grafico_crecimiento_var4 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var4, y = tasa_crecimiento_var4), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var4, y = tasa_tendencia_var4), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable4 IPIR_PAPEL: 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_var4)

Tras observar las tasas de crecimiento originales y tendencia se observa que, las tres variables tienen una tendencia muy marcada, sin embargo, es evidente la relación directa que tienen los siguientes indicadores: 1. ENER_V con IPIR_PAPEL: Se observa que la producción de cartón tiene un comportamiento similar al del consumo de energía de la región, esto indica que el sector del cartón está creciendo en producción en los mismos márgenes que las demás compañías manufactureras del valle del cauca, dado que, a mayor cantidad de energía requerida en los procesos, mayor es la producción. 2. Así mismo, el COR_V crece en las mismas proporciones del indicador anterior, y esto es explicado por que, a mayor producción de papel y cartón, mayor es el insumo disponible para la fabricación de empaques, siendo esta la materia prima clave para la fabricación de empaques.

El indicador que, si tiene un comportamiento un poco mas desajustado a lo anteriormente expuesto, pero no de manera negativa, sino positiva es el de las exportaciones de papel X_PAPEL, dado que, se observa que, aunque ha tenido volatibilidad entre 2019 y 2024, su tendencia va más al alza, lo cual, significa que gran parte de la producción nacional, se está exportando, siendo los principales mercados los centroamericanos y suramericanos.
Tras analizar el comportamiento de las gráficas originales vs estacionalidad, se observa que ninguno de las 4 tendencias presenta un ajuste representativo o significativo de estacionalidad sobre la original, esto se debe a que el componente estacional no es un factor que impacte o genere un dominio sobre ellas.

Modelo ARIMA

Modelo ARIMA para la variable 1,ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle

# 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

Identificación automática del modelo ARIMA

library(forecast)

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

# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)

## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.3214  -0.8831
## s.e.  0.1210   0.0789
## 
## sigma^2 = 145.5:  log likelihood = -593.63
## AIC=1193.27   AICc=1193.43   BIC=1202.34
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 1.00956 11.94281 8.333531 0.1633755 3.773785 0.7959477 -0.02951187

Estimación del modelo identificado automatico y validación de 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  0.321429   0.121041   2.6555  0.007918 ** 
## ma1 -0.883116   0.078896 -11.1934 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Ajuste del modelo ARIMA(1,1,1) 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, 1, 1))  # Especificamos directamente (p=1, d=1, q=1)  Ya que fue el resultado obtenido en el punto anterior.

# Mostrar resumen del modelo ajustado
summary(darima_auto)

## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.3214  -0.8831
## s.e.  0.1210   0.0789
## 
## sigma^2 = 145.5:  log likelihood = -593.63
## AIC=1193.27   AICc=1193.43   BIC=1202.34
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 1.00956 11.94281 8.333531 0.1633755 3.773785 0.7959477 -0.02951187

Validación de residuales o errores del modelo

# 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,1,1)
## Q* = 35.575, df = 22, p-value = 0.03371
## 
## Model df: 2.   Total lags used: 24

Pronóstico modelo ARIMA automático dentro de muestra o en el set de prueba

# 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

Interpretación modelo automatico (1,1,1):El gráfico muestra una clara diferencia entre lo pronosticado (azul) y lo real (rojo) frente a la variable ENER_V Consumo de energía de las grandes empresas industriales y comerciales (mercado no regulado) del Valle Pronóstico (Línea Azul): El modelo pronostica una tendencia gradual a la baja lo cual, representa una disminución del consumo de energía de empresas industriales del sector, lo cual, representa en nuestro sector de análisis una caída pequeña pero constante en la producción. Observado (Línea Roja): La realidad fue diferente, pues el consumo de energía cae aceleradamente, lo que indica una disminución abrupta en las operaciones de la cual aun no se recupera. En resumen, mientras el modelo preveía un proceso de decrecimiento paulatino, la operación real se vió afectada de una manera mucho mas agresiva. Esto refuerza la necesidad de mejorar el modelo de pronóstico para que considere otros factores que pueden afectar la producción de cartón.

Pronóstico automático dentro del set de prueba como tabla

# 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    238.68     234.3665
## 2 2024.833    232.39     233.9218
## 3 2024.917    228.59     233.7788

Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2025

Es decir, le sumamos al periodo de prueba una 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   234.3665
## 2 2024.833   233.9218
## 3 2024.917   233.7788
## 4 2025.000   233.7329

# 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 = 233.732900536782"

Modelo SARIMA para la variable 1:ENER_V

Identificación dautomática del modelo SARIMA

# Identificación automática modelo SARIMA
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,1)[12] 
## 
## Coefficients:
##           ma1    sma1
##       -0.7062  0.3214
## s.e.   0.0734  0.0801
## 
## sigma^2 = 136.9:  log likelihood = -589.52
## AIC=1185.04   AICc=1185.2   BIC=1194.11

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(0, 1, 1),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(0, 0, 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(0,1,1)(0,0,1)[12] 
## 
## Coefficients:
##           ma1    sma1
##       -0.7062  0.3214
## s.e.   0.0734  0.0801
## 
## sigma^2 = 136.9:  log likelihood = -589.52
## AIC=1185.04   AICc=1185.2   BIC=1194.11
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.5441839 11.58438 7.706434 -0.01019113 3.497339 0.7360527
##                    ACF1
## Training set 0.06840421

Validación de residuales del modelo automatico SARIMA

# 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,1)[12]
## Q* = 29.94, df = 22, p-value = 0.1199
## 
## Model df: 2.   Total lags used: 24

Pronóstico con el modelo SARIMA dentro del set de prueba-Gráfico líneas

# 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

Interpretación El resultado obtenido con este modelo si bien esta un poco mas ajustado a la realidad, se observa que al final suaviza un poco la caída, y de hecho podría verse una leve mejoría y esto obedece a que este tipo de modelos, combina información con rezagos y leves errores pasados, lo cual realiza que las caídas no sean tan bruscas y promedien con la tendencia previa y genere pronósticos más moderados y/o rebote. Esto toda vez que, al analizar el comportamiento similar en el pasado, es muy probable que se observaran leves mejorías de forma posterior a una crisis, por tanto, no prevé cambios estructurales. En este orden de ideas, se puede deducir que, el sector del papel y cartón va a guardar una tendencia al consumo de energía del sector manufacturero de la región, y pese a una caída en su crecimiento paulatina durante los últimos años, se prevé un rebote o mejora.

Pronóstico del modelo automático SARIMA en el set de prueba-Tabla

# 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    238.68     233.7317
## 2 2024.833    232.39     231.7294
## 3 2024.917    228.59     229.5779

Pronóstico del modelo automático SARIMA fuera de muestra, es decir, en enero 2025

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   233.7317
## 2 2024.833   231.7294
## 3 2024.917   229.5779
## 4 2025.000   234.1259

# 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 = 234.125893282209"

Interpretación El modelo ARIMA pronostica el uso de 233.732900536782 Gwh en el Valle del cauca para enero de 2025 lo cual, representa una mejora del 2,25% frente al mes inmediatamente anterior y del 0,48% frente al mismo periodo del año anterior. El modelo SARIMA pronostica el uso de 234.125893282209 Geh en el valle del cauca para enero de 2025, lo cual, representa una mejora del 2,5% frente al mes anterior y del 0,65% frente al año anterior. El pronostico mas alto es el SARIMA y el mas cercano al comportamiento de los periodos anteriores es el ARIMA, por lo cual, se concluye que el modelo ARIMA presenta un pronostico mucho mas conservador para el sector manufacturero.

Modelo ARIMA

Modelo ARIMA para la variable 2,COR_V Indicador de Exportaciones de papel y cartón

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

Identificación automática del modelo ARIMA

library(forecast)

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

# Mostrar el modelo seleccionado
summary(auto_arima_var2_no_seasonal)

## Series: train_ts_var2 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.3756  -0.9151
## s.e.  0.0899   0.0400
## 
## sigma^2 = 125.3:  log likelihood = -582.36
## AIC=1170.72   AICc=1170.88   BIC=1179.79
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
## Training set 1.001098 11.08276 8.390742 0.1434114 6.30082 0.8642991 0.01689634

Estimación del modelo identificado automatico y validación de Significancia de coeficientes

library(lmtest)
library(forcats)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_var2_no_seasonal)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.375603   0.089922   4.177 2.954e-05 ***
## ma1 -0.915110   0.039984 -22.887 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Ajuste del modelo ARIMA(1,1,1) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_var2_auto <- Arima(train_ts_var2, 
                order = c(1, 1, 1))  # Especificamos directamente (p=1, d=1, q=1)  Ya que fue el resultado obtenido en el punto anterior.

# Mostrar resumen del modelo ajustado
summary(darima_var2_auto)

## Series: train_ts_var2 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.3756  -0.9151
## s.e.  0.0899   0.0400
## 
## sigma^2 = 125.3:  log likelihood = -582.36
## AIC=1170.72   AICc=1170.88   BIC=1179.79
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
## Training set 1.001098 11.08276 8.390742 0.1434114 6.30082 0.8642991 0.01689634

Validación de residuales o errores del modelo

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_var2_auto)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 88.525, df = 22, p-value = 6.102e-10
## 
## Model df: 2.   Total lags used: 24

Pronóstico modelo ARIMA automático dentro de muestra o en el set de prueba

# Generar pronóstico para el conjunto de prueba
forecast_arima_var2 <- forecast(darima_var2_auto, h = length(test_ts_var2))  # Predecir los valores futuros

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

# Graficar pronóstico junto con los valores observados reales
p4_var2<- ggplot(forecast_data_var2, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado variable 2 COR_V") +
  xlab("Tiempo") + ylab("variable2: COR_v")

ggplotly(p4_var2)  # Convertir el gráfico en interactivo

Interpretación modelo automatico (1,1,1):

Tal y como ocurrió para la variable 1, el modelo Arima presenta una graficación estable y no refleja las caidas abruptas de los valores observados. Pronóstico automático dentro del set de prueba como tabla

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

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_var2 <- forecast(auto_arima_var2_no_seasonal, h = length(test_ts_var2))

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

# Mostrar la tabla
print(forecast_table_var2)

##     Tiempo Observado Pronosticado
## 1 2024.750     146.2     135.6101
## 2 2024.833     141.8     136.1398
## 3 2024.917     121.3     136.3387

Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2025

Es decir, le sumamos al periodo de prueba una 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_var2 <- forecast(auto_arima_var2_no_seasonal, h = length(test_ts_var2) + 1)

# Extraer el pronóstico del próximo trimestre
next_month_forecast_var2 <- data.frame(
  Tiempo = time(next_forecast_var2$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast_var2$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast_var2)

##     Tiempo Pronostico
## 1 2024.750   135.6101
## 2 2024.833   136.1398
## 3 2024.917   136.3387
## 4 2025.000   136.4134

# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_var2, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))

## [1] "Pronóstico para enero 2025: 2025 = 136.413430814363"

Modelo SARIMA para la variable 2:COR_

Identificación automática del modelo SARIMA

# Identificación automática modelo SARIMA
auto_arima_var2_no_seasonal <- auto.arima(train_ts_var2)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_var2_no_seasonal)

## Series: train_ts_var2 
## ARIMA(0,1,2)(1,0,1)[12] 
## 
## Coefficients:
##           ma1      ma2    sar1     sma1
##       -0.5804  -0.1917  0.8317  -0.4821
## s.e.   0.0835   0.0950  0.0987   0.1685
## 
## sigma^2 = 90.88:  log likelihood = -559.24
## AIC=1128.48   AICc=1128.89   BIC=1143.6

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima_var2 <- Arima(train_ts_var2, 
                order = c(0, 1, 2),  # (p,d,q) -> (1,0,1)
                seasonal = list(order = c(0, 0, 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(0,1,1)(0,0,1)[12] 
## 
## Coefficients:
##           ma1    sma1
##       -0.7062  0.3214
## s.e.   0.0734  0.0801
## 
## sigma^2 = 136.9:  log likelihood = -589.52
## AIC=1185.04   AICc=1185.2   BIC=1194.11
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.5441839 11.58438 7.706434 -0.01019113 3.497339 0.7360527
##                    ACF1
## Training set 0.06840421

Validación de residuales del modelo automatico SARIMA

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_var2)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,0,1)[12]
## Q* = 42.339, df = 21, p-value = 0.003814
## 
## Model df: 3.   Total lags used: 24

Pronóstico con el modelo SARIMA dentro del set de prueba-Gráfico líneas

# Generar pronóstico para el conjunto de prueba
forecast_arima_var2 <- forecast(darima_var2, h = length(test_ts_var2))  # Predecir los valores futuros

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

# 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

Pronóstico del modelo automático SARIMA en el set de prueba-Tabla

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

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_var2 <- forecast(auto_arima_var2_no_seasonal, h = length(test_ts_var2))

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

# Mostrar la tabla
print(forecast_table_var2)

##     Tiempo Observado Pronosticado
## 1 2024.750     146.2     136.7295
## 2 2024.833     141.8     137.4347
## 3 2024.917     121.3     125.4247

Pronóstico del modelo automático SARIMA fuera de muestra, es decir, en enero 2025

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_var2 <- forecast(auto_arima_var2_no_seasonal, h = length(test_ts_var2) + 1)

# Extraer el pronóstico del próximo mes
next_month_forecast_var2 <- data.frame(
  Tiempo = time(next_forecast_var2$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast_var2$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast_var2)

##     Tiempo Pronostico
## 1 2024.750   136.7295
## 2 2024.833   137.4347
## 3 2024.917   125.4247
## 4 2025.000   130.9372

# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month_var2 <- tail(next_month_forecast_var2, 1)
print(paste("Pronóstico para enero 2025:", next_month_var2$Tiempo, "=", next_month_var2$Pronostico))

## [1] "Pronóstico para enero 2025: 2025 = 130.937234177771"

Interpretación El modelo ARIMA pronostica un índice de COR_V de 136,41 para el valle del cauca, lo cual, está 6% por debajo de lo obtenido en enero de 2024, es decir, frente al mismo mes del año anterior. El modelo SARIMA pronostica un índice de COR_V de 130,937 para el valle del cauca, lo cual, está 10% por debajo de lo obtenido en enero de 2024, es decir, frente al mismo mes del año anterior. El pronostico mas alto es el ARIMA y es el mas cercano al comportamiento de los periodos anteriores por tanto, se concluye que el modelo ARIMA se ajusta mejor al comportamiento de años anteriores del sector del cartón.

Conclusión:

Dado que el pronostico es a la baja paulatina, se recomienda a la administración realizar uso eficiente de sus recursos físicos y así reducir el costo fijo de energía. Realizar planes de estructuración de costos, con la finalidad de generar eficiencias al interior de la compañía, invetir en procesos de automatización, con la finalidad de aumentar la producción con la menor cantidad de recursos. Si el mercado sigue a la baja, se recomienda ajustar la producción de conformidad a la demanda. Como estrategia comercial, se observa en la variable de exportaciones un crecimiento moderado, por lo cual, se recomienda realizar una importante gestión comercial en el exterior, crear nuevos productos y conseguir clientes con altos volúmenes de pedidos y hacer uso de incentivos dados por el gobierno nacional.

Extracción de señales y pronóstico ARIMA: Análisis financiero de los indicadores clave del sector de cartón y papel en Colombia.

Jéssica Torres

2025-08-31

Introducción

Metodología

Modelo ARIMA

Modelo SARIMA para la variable 1:ENER_V

Modelo ARIMA

Modelo SARIMA para la variable 2:COR_