Analisis Sector Externo del Valle
VANESSA
2025-04-02
#Cargar librerías necesarias
library(readxl) # Para leer archivos Excel
library(tseries) # Para pruebas de estacionariedad
library(forecast) # Para modelado ARIMA y pronósticos
library(ggplot2) # Para visualización de datos
library(plotly) # Para gráficos interactivos
library(timetk) library(readxl)
data_col <- read_excel("Base Taller2.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"))
View(data_col)Variable 1: TRM
# Convertir/declarar variable 1=TRM en serie de tiempo mensual
variable1_ts <- ts(data_col$TRM, start = c(2012, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 1847.516 1786.543 1767.793 1773.893 1795.907 1788.877 1784.903 1807.772
## 2013 1770.024 1794.198 1813.454 1830.091 1853.722 1910.188 1900.163 1904.411
## 2014 1964.111 2042.631 2017.022 1936.680 1915.203 1887.066 1858.366 1899.140
## 2015 2400.168 2422.994 2593.465 2492.766 2444.315 2561.667 2743.140 3027.064
## 2016 3287.158 3356.246 3124.169 2992.575 3002.174 2985.018 2971.474 2956.525
## 2017 2941.918 2880.985 2941.717 2876.237 2920.387 2966.522 3034.359 2971.017
## 2018 2866.496 2860.524 2842.387 2765.226 2863.577 2893.081 2885.532 2965.266
## 2019 3160.838 3112.360 3131.566 3158.608 3311.943 3250.837 3212.277 3419.289
## 2020 3321.721 3415.273 3908.308 3973.602 3853.472 3701.386 3659.323 3788.981
## 2021 3494.952 3551.997 3618.259 3654.614 3737.572 3687.401 3835.867 3882.250
## 2022 4002.746 3934.644 3799.549 3800.747 4015.570 3950.261 4389.145 4324.015
## 2023 4715.023 4822.552 4759.074 4530.264 4536.510 4205.262 4054.683 4077.895
## 2024 3922.445 3931.981 3899.277 3871.303 3867.880 4059.680 4038.110 4068.280
## Sep Oct Nov Dec
## 2012 1802.369 1805.740 1819.583 1791.352
## 2013 1917.924 1884.893 1923.471 1933.276
## 2014 1976.424 2048.908 2133.388 2351.757
## 2015 3074.955 2923.116 3007.221 3243.241
## 2016 2919.149 2938.747 3109.397 3003.960
## 2017 2917.878 2960.133 3014.114 2990.191
## 2018 3033.467 3092.217 3198.373 3216.620
## 2019 3403.978 3431.769 3412.555 3366.764
## 2020 3756.742 3832.831 3674.118 3461.582
## 2021 3823.625 3768.999 3906.607 3969.679
## 2022 4445.920 4726.851 4920.763 4793.103
## 2023 4007.122 4227.934 4035.887 3948.097
## 2024 4193.890 4266.350 4412.060 4386.390Variable 2: Exportaciones Valle
# Convertir/declarar variable 2=X_V en serie de tiempo mensual
variable2_ts <- ts(data_col$X_V, start = c(2012, 1), frequency = 12)
variable2_ts ## Jan Feb Mar Apr May Jun Jul
## 2012 159599364 180240665 191373990 161232070 198790367 172481442 194364495
## 2013 164853457 169801437 164553704 198509588 175509933 152859086 173179354
## 2014 120373398 156482020 160046060 184237919 178427402 152174190 192817659
## 2015 120882557 176969151 149977104 177977994 162833237 146026274 161483187
## 2016 128740720 189942224 167321121 147770050 150401389 124784925 94999973
## 2017 128273448 148968714 149447316 148640362 145481369 119980457 163823722
## 2018 132657531 160742577 147822277 172354823 143346837 132200919 142668441
## 2019 140010461 136674121 131806958 161263708 149089961 132475371 155450237
## 2020 139829192 135872501 136076573 120039763 135953465 135353513 163518257
## 2021 124798823 138082034 151649690 135855470 49767731 84260164 162444596
## 2022 108843084 144701812 157113802 141342682 130867409 152993170 134477390
## 2023 112288398 139351709 174700489 145292525 176664402 149265232 159393854
## 2024 117345624 140222685 141710961 157220167 161085481 142971963 149111281
## Aug Sep Oct Nov Dec
## 2012 202149989 192678930 217124571 182018485 187917021
## 2013 159632675 165283964 186448482 176526791 177119851
## 2014 210298339 200647922 195984020 203404691 172781130
## 2015 143577467 157897446 160756671 164555028 137107036
## 2016 155183649 195113834 222055829 190173966 201443946
## 2017 151509921 170903124 151608899 162331257 148773197
## 2018 168154486 139531023 159544847 163467561 140322064
## 2019 165578025 154506012 170845927 150933116 140487586
## 2020 153355466 169110909 158758029 155184516 171559349
## 2021 150574146 142640911 130498956 163384065 168378320
## 2022 144420197 137914582 141911619 161517130 162774116
## 2023 167957125 155852975 168219469 150665774 142357816
## 2024 148585302 154743385 180276265 153499206 172084853Variable 3: Importaciones Valle
# Convertir/declarar variable 1=M_V en serie de tiempo mensual
variable3_ts <- ts(data_col$M_V, start = c(2012, 1), frequency = 12)
variable3_ts## Jan Feb Mar Apr May Jun Jul
## 2012 355074079 556359400 408991672 344847316 411510051 390090495 394875972
## 2013 410666973 359139387 355342304 427014511 411642940 379639779 380236995
## 2014 412087737 374326560 413940374 471041989 470730946 421020081 496540003
## 2015 438172881 396410265 372288544 370972511 359980761 355598560 396045401
## 2016 327720078 308905727 317695275 345953466 300974616 289314046 282372332
## 2017 331357913 329140297 370590410 321930508 292526140 331623790 284971846
## 2018 333349548 323432567 313055506 347877391 335567028 329076039 350568656
## 2019 368632868 307594284 324104730 341831414 341779950 317345438 331304185
## 2020 321480804 386684901 299095593 304821736 286658639 268102810 326536395
## 2021 319889082 346588362 408436988 369834166 248993701 373647208 388686740
## 2022 435008827 393865542 486735774 424340748 421078617 420340822 320088971
## 2023 351172059 303315141 375266946 345832264 339176071 299875198 253045708
## 2024 352871095 330799782 320792544 572115119 363635093 307749248 326743464
## Aug Sep Oct Nov Dec
## 2012 427472544 366472004 421252726 413280397 356850071
## 2013 383357976 396772717 430514316 356858327 391124496
## 2014 455969188 460384545 533928375 431048217 417265246
## 2015 340035359 347999399 360625240 340066016 334401997
## 2016 378108884 326741403 324976024 350343625 321872761
## 2017 357719869 317298188 333985758 345585114 300749940
## 2018 364386573 271741552 381432147 341869649 303083199
## 2019 401498247 334206010 371487077 316632940 321711969
## 2020 297092358 296573918 332525591 349123453 376900177
## 2021 374419409 410794826 437442327 427677263 500470221
## 2022 479405715 453044681 419507074 360037022 399598223
## 2023 314414115 298036781 323767865 309797740 333716565
## 2024 357513502 335232735 366814925 412559322 366946648EXTRACCION DE SEÑALES
Extraccion Señales Variable 1: TRM
# 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: TRM",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extraccion de Señales Variable 2: Exportaciones 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: Exportaciones Valle",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extraccion de Señales 3: Importaciones 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: Importaciones Valle",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)TRM ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]Exportaciones Valle ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]Imporaciones Valle ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]Gráfico serie original VS ajustada Variable 1: TRM
# 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 = "orange", 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 TRM: 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: Exportaciones 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 = "orange", 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 X Valle: 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: Importaciones Valle
# 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 = "orange", 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 M Valle: 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)Grafico Serie Original VS Tendencia
*Tendencia Varibale 1: TRM
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" = "orange")) +
ggtitle("Variable 1 TRM: 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 Varibale 2: Exportaciones 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" = "orange")) +
ggtitle("Variable 2 X Valle: 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 Varibale 3: Imporatciones Valle
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" = "orange")) +
ggtitle("Variable 3 M Valle: 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)Calculo tasa de crecimiento de la serie de tendencia y original para la variable 1: TRM
#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] 144Gráfico variable original y tendencia variable 1 TRM: 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("Variable 1 TRM: 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)Gráfico variable original y tendencia variable 2 X valle: tasa de crecimiento anual
#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("Variable 2 X Valle: 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)#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("Variable 3 M Valle: 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)MODELO ARIMA
# Esta división idealmente podria se 80%-70% de los datos para entrenamiento y 20%-30% para prueba o test
# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2024 y el conjunto de prueba o test: noviembre 2024-diciembre 2024
train_size <- length(variable2_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable2_ts, end = c(2024, 9)) # Entrenamiento hasta septiembre 2024
test_ts <- window(variable2_ts, start = c(2024, 10)) # Prueba inicia desde oct2024Paso 1: Identificación del modelo
library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -5.0177, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 1 , una sola vez:
train_diff2 <- diff(train_ts, differences = 2) Diferenciación en niveles variable 2
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable2 = as.numeric(train_ts)), aes(x = Tiempo, y = variable2)) +
geom_line(color = "blue") +
ggtitle("Variable 2 Exportaciones Valle: 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)[-c(1,2)], variable2_Diff2 = as.numeric(train_diff2)), aes(x = Tiempo, y = variable2_Diff2)) +
geom_line(color = "red") +
ggtitle("variable 2 Exportaciones Valle: Serie Estacionaria (Dos diferenciación)") +
xlab("Tiempo") + ylab("variable 2 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivo# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_test_diff <- adf.test(train_diff2)
print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff2
## Dickey-Fuller = -10.715, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryIdentificación manual de p y q
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff2, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff2, 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 p=3 p=4 p=6 q=1* q6 (P optimo=1) (q óptimo=1)
El modelo óptimo para esta variable seria (3,2,6)
Paso 2. Estimación manual del modelo
Estimación del modelo identificado (3,2,6)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(3,2,6)) #Se va a estimar un modelo Arima de orden (3,2,6)
summary(manual_arima_model)## Series: train_ts
## ARIMA(3,2,6)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 ma4 ma5
## -0.8288 0.3319 0.7182 -0.7084 -1.4819 0.1354 1.6210 -0.1379
## s.e. 0.1096 0.1639 0.0983 0.0980 0.1356 0.1061 0.1053 0.1356
## ma6
## -0.4249
## s.e. 0.0993
##
## sigma^2 = 4.044e+14: log likelihood = -2757.37
## AIC=5534.75 AICc=5536.32 BIC=5564.92
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 828454 19372882 15084906 -1.31323 10.57041 0.7896183 -0.02725593Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.828775 0.109643 -7.5589 4.065e-14 ***
## ar2 0.331890 0.163852 2.0256 0.04281 *
## ar3 0.718156 0.098287 7.3067 2.738e-13 ***
## ma1 -0.708427 0.097972 -7.2309 4.798e-13 ***
## ma2 -1.481896 0.135614 -10.9273 < 2.2e-16 ***
## ma3 0.135362 0.106142 1.2753 0.20221
## ma4 1.620954 0.105306 15.3928 < 2.2e-16 ***
## ma5 -0.137938 0.135587 -1.0173 0.30899
## ma6 -0.424915 0.099346 -4.2771 1.893e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Paso 3. Validación de residuos del modelo estimado manual
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,2,6)
## Q* = 23.196, df = 15, p-value = 0.08007
##
## Model df: 9. Total lags used: 24Paso 4. Pronóstico (modelo manual)
Pronóstico en el test de prueba (oct, nov y dic 2024) y gráfico
#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva 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("Variable 2 Exportaciones Valle: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: 16826970cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 17256485A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic2024
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_manual <- forecast(manual_arima_model, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_manual <- data.frame(
Tiempo = time(arima_forecast_manual$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_manual$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_manual)## Tiempo Observado Pronosticado
## 1 2024.750 180276265 160849269
## 2 2024.833 153499206 142081870
## 3 2024.917 172084853 152448277Ahora 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 160849269
## 2 2024.833 142081870
## 3 2024.917 152448277
## 4 2025.000 150085238# 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 = 150085238.098813"