Series de tiempo: Extracción de señales y Modelo ARIMA
Geovani Zamora
2025-06-13
INTRODUCCIÓN
En este informe se hace un analisis de series temporales y pronostico de la empresa Bucanero S.A, una de las empresas con mayor participación en producción y ventas de Aves tanto a nivel nacional, como a nivel Valle del Cauca.
Pollos Bucaneros es una destaca empresa Colombiana, especializada en la producción y comercialización de productos alimenticios nutritivos, principalmente pollo. Es una de las avicolas lideres en Colombia y forma parte desde el año 2017 del gigante global de la industria Cargill.
Pollos Bucanero ofrece una alta gama de productos de productos de pollo, entre los que se encuentran, pollo fresco y congelado, producotos de valor agregado, opciones saborizadas y visceras
Instalar/Cargar librerias necesarias para el análisis
#Cargar librerías necesarias
library(readxl) # Para leer archivos Excel
library(tseries) # Para pruebas de estacionariedad
library(forecast) # Para modelado ARIMA y pronósticos
library(ggplot2) # Para visualización de datos
library(plotly) # Para gráficos interactivos
library(timetk) #timetk simplifica y acelera el análisis exploratorio, visualización, y preparación de datos temporales para modelado. Es ideal para quienes trabajan con series temporales en un flujo de trabajo "tidy" y buscan integrar análisis visuales, detección de patrones y forecasting en un solo paquete.Cargar base de datos
library(readxl)
data_col <- read_excel("C:/Users/Geovani Zamora/OneDrive/Escritorio/CASO 2 R/Datos Caso 2.xlsx", col_types = c("date",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric"))JUSTIFICACIÓN
Las variables seleccionadas para efectur el analisis y pronostico de los desempeño de la compañia son Variable 1 importación de cereales, Variable 2 producción de Pollos a nivel nacional y Varible 3 (producción de Pollos en el Valle del Cauca).
Las tres variables seleccionadas en primera instancia, se podria indicar que presentar una relacion directa y fuerte, con una dependencia considerable. Dado, que segun cifras de la ANDI, la industria avicola es el mayo consumidor de alimentos para alimentos, representando el 65% del total de importaciones de alimentos para animales, destacando el consumo de maiz amarillo, cereal con mayor uso en al alimentación de las aves.
PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):
Variable Importación de Cereales
# Convertir/declarar variable 1=Importación de Cereales en serie de tiempo mensual
variable1_ts <- ts(data_col$M_CEREAL, start = c(2012, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul
## 2012 130488.68 170760.24 150305.77 112970.29 171028.68 173377.41 130779.12
## 2013 174878.60 148872.59 102558.49 203834.71 143799.75 140773.61 228355.85
## 2014 105236.84 125248.74 146145.21 238881.10 231599.26 179125.59 92573.43
## 2015 181346.01 138858.00 166598.35 180075.91 199772.41 98064.31 96163.32
## 2016 140790.90 136539.24 153555.26 245910.28 157028.00 57305.22 88277.80
## 2017 145475.11 116937.84 214756.32 204681.76 75436.90 97989.00 96763.43
## 2018 154842.08 99009.76 150132.47 188181.24 182990.25 156560.18 104369.47
## 2019 182837.33 122823.96 143154.85 187626.02 181138.52 114712.02 158222.13
## 2020 226081.46 206790.11 133135.30 194238.41 195012.50 140688.95 127811.57
## 2021 168362.03 175159.63 241602.31 216342.77 203543.60 214196.52 189115.07
## 2022 231997.90 176420.62 377815.90 275554.59 316034.78 330913.11 267161.17
## 2023 232492.41 253829.89 336580.50 248993.64 220930.55 185685.56 206086.04
## 2024 202870.30 179162.77 217248.34 236529.82 262284.36 177963.37 123368.40
## Aug Sep Oct Nov Dec
## 2012 209066.36 154539.70 162266.56 177411.83 179662.26
## 2013 165477.74 129904.20 146619.12 136963.24 130192.23
## 2014 106245.69 155162.32 145195.89 75898.77 88385.75
## 2015 121567.49 133577.01 114400.86 108865.64 110490.03
## 2016 119695.30 136907.56 133528.75 100765.46 84857.12
## 2017 112388.61 140841.41 125074.25 92688.23 87718.28
## 2018 137571.35 113941.63 145878.28 124574.62 109430.60
## 2019 239728.89 168349.34 114161.10 146884.46 116240.81
## 2020 166267.17 133404.40 144408.08 163827.69 130209.92
## 2021 236088.30 191084.91 303131.61 203863.59 227212.86
## 2022 317699.79 242911.52 273134.00 241116.53 273212.27
## 2023 283608.77 161500.03 183800.39 163547.31 252087.78
## 2024 234965.51 166055.79 208222.38 199493.94 205484.87Variable Producción de Pollos a Nivel Nacional
# Convertir/declarar la Variable Pollo en serie de tiempo mensual
variable2_ts <- ts(data_col$POLLO, start = c(2012, 1), frequency = 12)
variable2_ts## Jan Feb Mar Apr May Jun Jul
## 2012 91721.61 94142.30 88748.18 92013.41 93279.24 91314.75 85934.76
## 2013 102363.98 102875.96 98330.23 96164.43 102075.96 108350.85 104104.91
## 2014 106197.36 110133.59 103197.17 105954.06 109403.32 112677.66 111061.74
## 2015 115792.62 118873.07 113713.78 119043.57 120466.53 113405.45 113300.36
## 2016 124206.78 120032.40 117024.65 121569.62 119589.96 120656.25 119707.67
## 2017 131659.80 130484.69 121952.89 123962.79 122845.66 126745.99 130163.07
## 2018 126092.16 128869.53 125704.17 133765.31 137389.12 141774.82 137681.85
## 2019 134114.63 135908.08 129794.01 131683.46 140610.05 147000.19 139211.72
## 2020 143717.95 147147.06 137739.46 136479.81 116092.75 102369.20 118464.78
## 2021 133650.94 134213.48 129669.19 145854.17 147160.96 128269.37 129431.15
## 2022 148780.18 147752.77 139974.80 148835.06 150382.43 152736.19 149930.18
## 2023 152778.62 154942.57 141342.54 146778.12 138725.93 143836.02 151160.53
## 2024 146670.25 151062.12 145678.30 146285.84 152650.08 157002.95 148202.52
## Aug Sep Oct Nov Dec
## 2012 90465.71 95904.98 91873.01 97307.48 99554.95
## 2013 109811.58 113429.54 110209.13 114715.85 111837.83
## 2014 116671.17 116616.51 121499.06 126515.35 119226.47
## 2015 118715.14 119359.22 119806.21 125709.60 126202.13
## 2016 118099.12 126959.91 130587.16 129761.49 130727.74
## 2017 128387.97 136740.14 134601.16 138998.94 137064.01
## 2018 138700.74 143771.77 134467.86 140107.07 141334.79
## 2019 139890.73 149043.48 143802.26 152612.96 149506.78
## 2020 138648.94 143973.52 144064.42 148785.46 142301.21
## 2021 141917.90 150464.44 149214.12 150818.57 153670.57
## 2022 149606.80 156882.57 160015.22 158998.13 156231.77
## 2023 152113.55 160168.91 156991.78 161527.41 156763.90
## 2024 152107.76 156103.22 151683.90 158293.69 156027.51Variable Producción de Pollos en el Valle del Cauca
# Convertir/declarar las ventas de POLLO en el Valle en serie de tiempo mensual
variable3_ts <- ts(data_col$POLLO_V, start = c(2012, 1), frequency = 12)
variable3_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 13518.16 13579.57 13502.11 14575.28 14565.48 13434.11 12662.66 13002.30
## 2013 14003.60 13598.02 17994.93 15756.23 15247.82 15302.08 15560.96 16787.54
## 2014 17160.85 18042.63 17530.63 18583.41 18531.03 19541.67 19007.82 19535.39
## 2015 21967.08 22318.20 20011.13 21414.71 21776.54 21331.25 20859.23 21857.07
## 2016 25419.39 23836.67 23425.54 25223.50 24874.07 24428.70 24432.43 24351.71
## 2017 26767.99 27349.59 25185.77 25243.04 25288.09 26999.44 27039.38 27198.80
## 2018 26888.51 27943.05 26847.86 28596.43 28700.71 32243.89 31172.15 30795.92
## 2019 28648.63 28152.79 27272.11 27203.26 28599.35 30283.57 29120.47 29407.24
## 2020 29349.13 29540.82 27833.77 25791.75 21732.97 18630.33 22950.05 27290.86
## 2021 25884.72 25811.36 24460.44 27693.47 27901.31 17622.96 20657.57 27045.41
## 2022 29200.71 30401.75 26620.46 28501.35 27980.21 30308.21 28577.47 30309.97
## 2023 30902.26 30673.37 26606.57 29190.77 26320.95 28955.32 29441.06 29927.85
## 2024 31190.50 30788.20 29047.01 28744.59 30447.01 30274.53 30194.80 30826.10
## Sep Oct Nov Dec
## 2012 13811.78 12623.49 14419.18 13997.19
## 2013 17165.77 16901.55 18339.54 17439.61
## 2014 18206.56 19826.65 22212.84 20839.27
## 2015 21643.92 22028.13 22155.39 25570.60
## 2016 26439.72 26252.85 26615.81 26619.36
## 2017 29612.29 29085.95 28658.94 29722.71
## 2018 31105.56 28728.88 29415.70 28665.26
## 2019 30889.64 28887.32 30416.51 29872.71
## 2020 28458.63 28632.59 29462.15 27626.71
## 2021 29695.43 29012.25 29864.17 30304.56
## 2022 30967.86 33589.17 31595.13 32182.50
## 2023 31185.45 31502.72 29734.15 31056.28
## 2024 32768.56 32987.00 33245.00 33140.00Extracción de señales
Gráfico inicial de la variable Importación de Cereales en niveles -Original
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("2012-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)Extracción señales variable Importación de Cereales
# 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)Extracción señales variable Producción de Pollo a Nivel Nacional
# 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",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extracción señales variable Producción de Pollo en el Valle
# 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 = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Se crea la variable1 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]
variable1_sa## Jan Feb Mar Apr May Jun Jul
## 2012 124003.77 182286.03 125191.30 72932.71 146333.62 185436.31 155366.06
## 2013 168393.70 160398.38 77444.03 163797.13 119104.69 152832.51 252942.79
## 2014 98751.93 136774.53 121030.75 198843.52 206904.19 191184.49 117160.37
## 2015 174861.11 150383.79 141483.89 140038.33 175077.35 110123.20 120750.26
## 2016 134306.00 148065.03 128440.79 205872.70 132332.94 69364.11 112864.74
## 2017 138990.21 128463.63 189641.86 164644.18 50741.84 110047.90 121350.37
## 2018 148357.18 110535.55 125018.00 148143.66 158295.19 168619.07 128956.41
## 2019 176352.43 134349.75 118040.38 147588.45 156443.46 126770.92 182809.07
## 2020 219596.55 218315.90 108020.83 154200.83 170317.44 152747.85 152398.51
## 2021 161877.13 186685.42 216487.85 176305.19 178848.54 226255.42 213702.01
## 2022 225512.99 187946.41 352701.43 235517.01 291339.72 342972.01 291748.11
## 2023 226007.51 265355.68 311466.04 208956.07 196235.48 197744.45 230672.98
## 2024 196385.40 190688.56 192133.88 196492.25 237589.30 190022.27 147955.34
## Aug Sep Oct Nov Dec
## 2012 192426.40 170789.35 165903.84 201941.46 200046.05
## 2013 148837.77 146153.85 150256.40 161492.87 150576.02
## 2014 89605.72 171411.97 148833.17 100428.40 108769.54
## 2015 104927.52 149826.65 118038.14 133395.27 130873.82
## 2016 103055.34 153157.21 137166.03 125295.09 105240.91
## 2017 95748.64 157091.06 128711.53 117217.86 108102.07
## 2018 120931.38 130191.28 149515.56 149104.25 129814.39
## 2019 223088.93 184598.99 117798.38 171414.09 136624.60
## 2020 149627.20 149654.04 148045.36 188357.32 150593.71
## 2021 219448.33 207334.55 306768.89 228393.22 247596.64
## 2022 301059.82 259161.17 276771.28 265646.15 293596.06
## 2023 266968.80 177749.67 187437.67 188076.94 272471.57
## 2024 218325.54 182305.44 211859.66 224023.57 225868.66Se crea la variable2 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]
variable2_sa## Jan Feb Mar Apr May Jun Jul
## 2012 92409.75 93814.55 95270.61 94628.05 96081.33 94765.75 90284.85
## 2013 103052.12 102548.21 104852.66 98779.07 104878.04 111801.86 108454.99
## 2014 106885.49 109805.84 109719.60 108568.71 112205.40 116128.66 115411.83
## 2015 116480.76 118545.32 120236.22 121658.22 123268.61 116856.45 117650.45
## 2016 124894.91 119704.64 123547.09 124184.26 122392.05 124107.26 124057.75
## 2017 132347.93 130156.93 128475.32 126577.43 125647.74 130197.00 134513.16
## 2018 126780.29 128541.77 132226.60 136379.96 140191.20 145225.83 142031.94
## 2019 134802.77 135580.33 136316.44 134298.10 143412.13 150451.19 143561.80
## 2020 144406.08 146819.31 144261.90 139094.45 118894.84 105820.21 122814.86
## 2021 134339.07 133885.72 136191.62 148468.81 149963.04 131720.37 133781.23
## 2022 149468.32 147425.02 146497.23 151449.71 153184.51 156187.19 154280.27
## 2023 153466.75 154614.81 147864.97 149392.76 141528.01 147287.02 155510.62
## 2024 147358.38 150734.37 152200.73 148900.48 155452.16 160453.96 152552.61
## Aug Sep Oct Nov Dec
## 2012 90733.65 90744.67 88656.81 90193.87 94676.51
## 2013 110079.51 108269.23 106992.93 107602.24 106959.39
## 2014 116939.11 111456.21 118282.85 119401.74 114348.03
## 2015 118983.08 114198.91 116590.01 118595.99 121323.69
## 2016 118367.06 121799.61 127370.96 122647.88 125849.30
## 2017 128655.91 131579.83 131384.95 131885.33 132185.56
## 2018 138968.68 138611.46 131251.66 132993.46 136456.34
## 2019 140158.66 143883.17 140586.06 145499.35 144628.34
## 2020 138916.87 138813.21 140848.21 141671.85 137422.77
## 2021 142185.84 145304.13 145997.91 143704.96 148792.12
## 2022 149874.74 151722.26 156799.01 151884.52 151353.32
## 2023 152381.49 155008.60 153775.58 154413.79 151885.45
## 2024 152375.70 150942.92 148467.69 151180.08 151149.06Se crea la variable3 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]
variable3_sa## Jan Feb Mar Apr May Jun Jul Aug
## 2012 13246.13 13339.90 14589.57 14989.95 15440.57 14607.12 13754.40 12909.68
## 2013 13731.57 13358.34 19082.40 16170.90 16122.90 16475.09 16652.71 16694.93
## 2014 16888.82 17802.95 18618.09 18998.08 19406.11 20714.68 20099.57 19442.77
## 2015 21695.04 22078.52 21098.59 21829.38 22651.62 22504.26 21950.98 21764.45
## 2016 25147.36 23597.00 24513.01 25638.17 25749.15 25601.71 25524.17 24259.10
## 2017 26495.96 27109.91 26273.24 25657.71 26163.18 28172.45 28131.13 27106.19
## 2018 26616.47 27703.37 27935.33 29011.10 29575.80 33416.90 32263.90 30703.31
## 2019 28376.59 27913.11 28359.57 27617.93 29474.44 31456.58 30212.22 29314.62
## 2020 29077.10 29301.15 28921.23 26206.42 22608.05 19803.34 24041.79 27198.25
## 2021 25612.68 25571.68 25547.90 28108.14 28776.39 18795.97 21749.32 26952.79
## 2022 28928.67 30162.08 27707.93 28916.02 28855.30 31481.22 29669.21 30217.36
## 2023 30630.22 30433.69 27694.04 29605.44 27196.03 30128.33 30532.80 29835.23
## 2024 30918.47 30548.52 30134.47 29159.26 31322.10 31447.54 31286.55 30733.48
## Sep Oct Nov Dec
## 2012 12769.04 11840.06 13282.08 12922.82
## 2013 16123.03 16118.11 17202.43 16365.24
## 2014 17163.82 19043.21 21075.74 19764.90
## 2015 20601.18 21244.69 21018.29 24496.23
## 2016 25396.98 25469.42 25478.70 25544.99
## 2017 28569.55 28302.51 27521.84 28648.35
## 2018 30062.82 27945.44 28278.60 27590.89
## 2019 29846.90 28103.89 29279.41 28798.34
## 2020 27415.90 27849.16 28325.04 26552.35
## 2021 28652.69 28228.81 28727.07 29230.19
## 2022 29925.12 32805.73 30458.03 31108.13
## 2023 30142.71 30719.29 28597.04 29981.91
## 2024 31725.83 32203.56 32107.90 32065.63Gráfico serie original VS ajustada Variable 1
# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var1 <- ggplot() +
geom_line(aes(x = fechas_var1, y = variable1_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("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)Gráfico serie original VS ajustada Variable 2
# Crear vector de fechas correctamente alineado con la serie
fechas_var2 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("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
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)Gráfico serie original VS ajustada Variable 3
# Crear vector de fechas correctamente alineado con la serie
fechas_var3 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var3 <- ggplot() +
geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("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
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)Ahora graficamos serie original vs tendencia
Primero se debe obtener la tendencia de cada variable y luego graficarla
Tendencia Variable 1
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("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)Tendencia Variable 2
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("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)Tendencia Variable 3
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("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)Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:
Tasa de crecimiento de la serie de tendencia y original para la variable 1
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var1 <- (variable1_ts[(13:length(variable1_ts))] / variable1_ts[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1 <- (tendencia_var1[(13:length(tendencia_var1))] / tendencia_var1[1:(length(tendencia_var1) - 12)] - 1) * 100
# Crear vector de fechas corregido, es decir que inicie desde enero 2013
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))
# Verificar longitudes
print(length(fechas_corregidas_var1))## [1] 144print(length(tasa_crecimiento_var1))## [1] 144print(length(tasa_tendencia_var1))## [1] 144*Gráfico variable original y tendencia variable 1: tasa de crecimiento anual**
library(ggplot2)
library(plotly)
# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("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)Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))
# Verificar longitudes
print(length(fechas_corregidas_var2))## [1] 144print(length(tasa_crecimiento_var2))## [1] 144print(length(tasa_tendencia_var2))## [1] 144# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("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)Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))
# Verificar longitudes
print(length(fechas_corregidas_var3))## [1] 144print(length(tasa_crecimiento_var3))## [1] 144print(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)Analizar la tasa de crecimiento anual ayuda a detectar cambios en el entorno económico que afectan el sector. Se pueden prever crisis o períodos de auge y prepararse para ellos.
Modelo ARIMA
División en conjunto de entrenamiento y prueba para la variable 3 (Producción de Pollo en el Valle) que es la elegida para pronosticar
El código siguiente divide una serie temporal (variable3_ts) en dos subconjuntos:
Conjunto de entrenamiento (train): Datos desde enero de 2012 hasta septiembre de 2024. Conjunto de prueba (test): Datos desde octubre de 2024 hasta diciembre de 2024.
Esto se hace para evaluar el desempeño de modelos de predicción en datos no vistos.
train_size <- length(variable3_ts) - 12 # Se deja fuera los últimos 12 valores para usarlos como set de prueba.
train_ts <- window(variable3_ts, end = c(2023, 12)) # Entrenamiento hasta noviembre 2024
test_ts <- window(variable3_ts, start = c(2024, 1)) # Prueba inicia desde diciembre 2024Modelo ARIMA automático normal (sin tener en cuenta el factor estacional)
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(2,1,2) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 drift
## 1.5763 -0.7390 -1.8570 0.9340 112.0795
## s.e. 0.0716 0.0699 0.0451 0.0471 66.5356
##
## sigma^2 = 2951659: log likelihood = -1266.5
## AIC=2545 AICc=2545.62 BIC=2562.78
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -2.288793 1681.866 1132.138 -0.288442 4.816199 0.4622445
## ACF1
## Training set -0.04520703Estimació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 1.576286 0.071616 22.0103 < 2e-16 ***
## ar2 -0.739024 0.069882 -10.5754 < 2e-16 ***
## ma1 -1.857007 0.045123 -41.1547 < 2e-16 ***
## ma2 0.933992 0.047054 19.8493 < 2e-16 ***
## drift 112.079516 66.535564 1.6845 0.09208 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El modelo muuestra una alta significancia
# Ajuste del modelo ARIMA(1,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(2, 1, 2)) # Especificamos directamente (p=4, d=1, q=2)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 1.5912 -0.7417 -1.8669 0.9444
## s.e. 0.0688 0.0697 0.0391 0.0401
##
## sigma^2 = 2983157: log likelihood = -1267.8
## AIC=2545.61 AICc=2546.04 BIC=2560.42
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 216.1118 1696.931 1158.031 0.6512176 4.887072 0.4728163
## ACF1
## Training set -0.04851695Validació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(2,1,2)
## Q* = 28.522, df = 20, p-value = 0.0976
##
## Model df: 4. Total lags used: 24Pronóstico modelo ARIMA automático dentro de muestra o en el set de prueba (2,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("variable3")
ggplotly(p4auto) # Convertir el gráfico en interactivoInterpretación modelo automatico (2,1,2):El modelo automático (2,1,2) parece pronosticar de forma aceptable dentro de la prueba, arrojando en la prueba pronostico que se encuentran por debajo de los observados. se percibe captura de los los puntos de quiebre. Es un modelo tentativo adecuado para pronpostico fuera de muestra o a futuro.
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.000 31190.50 30225.75
## 2 2024.083 30788.20 29534.82
## 3 2024.167 29047.01 29077.75
## 4 2024.250 28744.59 28886.12
## 5 2024.333 30447.01 28940.09
## 6 2024.417 30274.53 29185.01
## 7 2024.500 30194.80 29549.44
## 8 2024.583 30826.10 29961.12
## 9 2024.667 32768.56 30358.96
## 10 2024.750 32987.00 30700.07
## 11 2024.833 33245.00 30961.98
## 12 2024.917 33140.00 31140.98Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2025
# 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.000 30225.75
## 2 2024.083 29534.82
## 3 2024.167 29077.75
## 4 2024.250 28886.12
## 5 2024.333 28940.09
## 6 2024.417 29185.01
## 7 2024.500 29549.44
## 8 2024.583 29961.12
## 9 2024.667 30358.96
## 10 2024.750 30700.07
## 11 2024.833 30961.98
## 12 2024.917 31140.98
## 13 2025.000 31247.82# 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 = 31247.8153108169"CONCLUSIÓN
El modelo automático ARIMA(2,1,2) en un pronostico para el mes de enero de 2025, presenta un valor superior al pronostico de prueba arrojado para el mes de diciemnbre de 2024, pero se encuentra por debajo del dato observado para el mes de diciembre de 2024. Destacó por su mayor precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.