1 Pronostico para la demanda de energia en Sincelejo

1.1 Punto A

library(readxl)
energy <- read_xlsx("energy.xlsx")
energy_dem <- energy %>% 
  summarise_by_time(.date_var = Fecha, .by = "month",
                                        demanda_real = mean(demanda_real, .na.rm = TRUE))
serie <- ts (energy_dem$demanda_real, start = c(2000, 1), end = c(2023, 9), 
              frequency = 12)
ts_plot(serie, slider = T,
        color = "blue",
        title = "Demanda mensual de energia en sincelejo 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")

la serie de energía en sincelejo sube con el tiempo, sobre todo después de 2020 donde hay un salto fuerte, la serie tiene un patrón que se repite cada año. Eso muestra que hay estacionalidad, los efectos se ven en meses específicos. En diciembre y enero la demanda sube por fiestas y vacaciones, a mitad de año también puede subir por vacaciones escolares. La estacionalidad pasa porque la energía depende del clima, fiestas y la actividad económica. cuando hace más calor la gente usa más ventiladores y aires. en diciembre y mitad de año se consume más en hogares y comercios.

ts_decompose(serie, type = "both")

La descomposición aditiva asume que la estacionalidad es igual en todo el periodo. la multiplicativa ajusta la estacionalidad al nivel de la serie. como la variación crece cuando la serie crece, la forma multiplicativa describe mejor el comportamiento.

1.2 Punto B

demanda_descomp <- decompose(serie, type = "additive")
demanda_descomp <- serie - demanda_descomp$seasonal
demanda <- data.frame("demanda_ob" = as.numeric(serie)) %>% 
 mutate(demanda_descomp = as.numeric(demanda_descomp)) %>%
  mutate(Fecha = fechas <- seq(
  from = as.Date("2000-01-01"),
  to   = as.Date("2023-09-01"),
  by   = "month"
))
ts_plot(demanda, slider = T,
      title = "Demanda mensual de energia en sincelejo normal y desestacionalizada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")

La serie desestacionalizada elimina los picos que se repetían en los mismos meses cada año. en la gráfica las dos series son casi iguales, pero la línea naranja (sin estacionalidad) es más suave. eso muestra que el efecto calendario se redujo bastante. antes se veían subidas en meses como diciembre o mitad de año. después de quitar la estacionalidad esas variaciones desaparecen y queda solo la tendencia y la parte aleatoria,con esto se entiende mejor el crecimiento real de la demanda sin la influencia de las fiestas, el clima o las vacaciones.

1.3 Punto C

1.3.1 Correlograma con tres años de rezago

acf(serie, 
    main = "Correlograma ACF - Demanda Energía", 
    lag.max = 36)

El correlograma muestra cómo la serie se correlaciona consigo misma en diferentes rezagos. En este caso, la autocorrelación disminuye de forma lenta a medida que aumentan los lags. Esto indica que los valores pasados siguen influyendo de manera prolongada en los valores futuros, lo cual es un signo de no estacionariedad.

1.3.2 Prueba Dickey Fuller

dickey_full <- ur.df(serie, type = "trend", lags = 0)
summary(dickey_full)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.699  -2.332  -0.471   1.914  53.296 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.349890   0.862926  -0.405   0.6854  
## z.lag.1     -0.016088   0.011675  -1.378   0.1693  
## tt           0.016660   0.007653   2.177   0.0303 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.97 on 281 degrees of freedom
## Multiple R-squared:  0.01712,    Adjusted R-squared:  0.01012 
## F-statistic: 2.447 on 2 and 281 DF,  p-value: 0.0884
## 
## 
## Value of test-statistic is: -1.3781 2.2545 2.4469 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36
dwtest(dickey_full@testreg)
## 
##  Durbin-Watson test
## 
## data:  dickey_full@testreg
## DW = 2.0209, p-value = 0.5228
## alternative hypothesis: true autocorrelation is greater than 0
par(mfrow = c(1, 2))
acf(serie, )
pacf(serie)

serie_dif <- diff(serie, differences = 1)

La prueba Dickey-Fuller contrasta la hipótesis de raíz unitaria. Al usar un rezago igual a cero, el resultado no rechaza la hipótesis nula de no estacionariedad. Esto confirma estadísticamente lo que ya sugería el correlograma: la serie no es estacionaria. Para volverla estacionaria sería necesario aplicar el metodo de diferenciación.

1.4 Punto D

dic_full_diff <- ur.df(serie_dif, type = "trend", lags = 0)
options(scipen = 999)
summary(dic_full_diff)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.168  -2.159  -0.369   1.873  54.277 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) -0.747230   0.835254  -0.895               0.372    
## z.lag.1     -1.020706   0.059707 -17.095 <0.0000000000000002 ***
## tt           0.009245   0.005119   1.806               0.072 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.999 on 280 degrees of freedom
## Multiple R-squared:  0.5107, Adjusted R-squared:  0.5072 
## F-statistic: 146.1 on 2 and 280 DF,  p-value: < 0.00000000000000022
## 
## 
## Value of test-statistic is: -17.0951 97.4168 146.1252 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36
par(mfcol = c(1, 2))
acf(serie_dif)
pacf(serie_dif)

Al diferenciar la serie una vez, se volvió estacionaria. eso significa que d = 1. Cuando revisamos el correlograma de la serie diferenciada, la acf corta rápido en el rezago 1 y luego se mantiene dentro de las bandas. eso es típico de un ma(1). la pacf no muestra rezagos significativos más allá de cero, por lo que no hay evidencia de un ar(p).

p = 0, d = 1, q = 1. el modelo sugerido es arima(0,1,1)

1.4.1 Estimar modelo arima (0,1,1)

  arima <- arima(serie, c(0, 1, 1), seasonal = list(order = c(0,1,1), perid = 12))
  stargazer(arima, type = "HTML")

1.5 Punto E

demanda <- demanda %>% mutate(resid = arima$residuals,
                              demanda_est = as.numeric(demanda_ob-resid)) %>% 
  select(demanda_ob, demanda_est, Fecha)
ts_plot(demanda, slider = T,
      title = "Demanda mensual de energia en sincelejo normal y ajustada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")

Muestra la demanda mensual de energía en sincelejo entre 2000 y 2023, comparando los valores observados con los estimados por el modelo arima(0,1,1). la serie real y la estimada se ajustan bastante bien, sobre todo después de 2010, lo que indica que el modelo capta la tendencia de crecimiento de la demanda.

par(mfcol = c(1, 2))
hist(arima$residuals)
qqnorm(arima$residuals)
qqline(arima$residuals, col = "blue")

jarque.bera.test(arima$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima$residuals
## X-squared = 6658.9, df = 2, p-value < 0.00000000000000022

Muestra el análisis de los residuos del arima(0,1,1). el histograma refleja que están centrados en cero, pero concentrados en un rango pequeño. el gráfico Normal Q-Q Plot muestra que los residuos no siguen bien la normalidad, ya que hay desviaciones en las colas. esto indica que el modelo no cumple el supuesto de normalidad en los errores.

arima_mod <- auto.arima(serie, trace = T)
## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)(1,0,1)[12] with drift         : Inf
##  ARIMA(0,1,0)            with drift         : 1910.797
##  ARIMA(1,1,0)(1,0,0)[12] with drift         : 1924.82
##  ARIMA(0,1,1)(0,0,1)[12] with drift         : 1913.399
##  ARIMA(0,1,0)                               : 1910.62
##  ARIMA(0,1,0)(1,0,0)[12] with drift         : 1921.764
##  ARIMA(0,1,0)(0,0,1)[12] with drift         : 1911.378
##  ARIMA(0,1,0)(1,0,1)[12] with drift         : 1918.505
##  ARIMA(1,1,0)            with drift         : 1913.599
##  ARIMA(0,1,1)            with drift         : 1912.805
##  ARIMA(1,1,1)            with drift         : 1915.305
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,0)                               : 1914.516
## 
##  Best model: ARIMA(0,1,0)
demanda <- demanda %>% mutate(resid = arima_mod$residuals,
                              demanda_est = as.numeric(demanda_ob-resid)) %>% 
  select(demanda_ob, demanda_est, Fecha)
ts_plot(demanda, slider = T,
      title = "Demanda mensual de energia en sincelejo normal y ajustada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")

Muestra la misma serie de demanda mensual, ahora ajustada con el modelo seleccionado por la función auto.arima. en este caso, la línea de la demanda estimada (demanda_descomp) se ajusta mejor a los valores observados, lo que sugiere que el nuevo modelo representa mejor el comportamiento de la serie.

par(mfcol = c(1, 2))
hist(arima_mod$residuals)
qqnorm(arima_mod$residuals)
qqline(arima_mod$residuals, col = "blue")

jarque.bera.test(arima_mod$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima_mod$residuals
## X-squared = 9670.1, df = 2, p-value < 0.00000000000000022
options(scipen = 99999)

Muestra los residuos del modelo auto.arima. aunque siguen centrados en cero, el Normal Q-Q Plot indica que todavía hay colas pesadas y que no hay normalidad perfecta. sin embargo, comparado con el primer modelo, los residuos están mejor ajustados y el modelo es más confiable para pronosticar.

1.5.1 Pronostico de los siguientes 5 meses

valores_fut <- forecast(arima_mod, h = 5, level = c(0.95))
print(valores_fut)
##          Point Forecast    Lo 95    Hi 95
## Oct 2023       229.9245 216.1741 243.6749
## Nov 2023       229.9245 210.4785 249.3705
## Dec 2023       229.9245 206.1081 253.7409
## Jan 2024       229.9245 202.4237 257.4253
## Feb 2024       229.9245 199.1776 260.6714
par(mfcol = c(1, 1))
plot(valores_fut)

Muestra el pronóstico de la demanda de energía para los 5 meses siguientes a la última observación. la línea azul representa los valores pronosticados y la franja gris el intervalo de confianza al 95%. los resultados muestran que la demanda se mantendrá estable 229 Mwh.

2 Pronostico para la demanda de pasajeros

2.1 Punto A

data("AirPassengers")
ts_plot(AirPassengers, slider = T,
        title = "Demanda del número de pasajeros de aerolíneas internacionales",
        Ytitle= "Número de pasajeros",
        Xtitle= "Año-mes")

Notamos que si se evidencia estacionalidad, en la gráfica se ve un patrón que se repite cada año y es que sube fuerte en los meses de verano (junio, julio, agosto) y baja en invierno (noviembre, diciembre, enero). Esto pasa porque la demanda de vuelos internacionales crece en vacaciones de verano y también en diciembre por navidad. en esos meses la gente viaja más por turismo y reuniones familiares. en cambio, en los meses fríos de inicio de año (enero-febrero) y después del verano (septiembre-octubre) la demanda baja porque hay menos vacaciones y menos viajes de ocio.

ts_decompose(AirPassengers, type = "both")

Aqui notamos una estacionalidad clara, en la parte de estacional aparecen picos y caídas que se repiten todos los años en los mismos meses.sube en verano (junio-agosto) y diciembre. Baja en los primeros meses del año (enero-febrero) y en otoño (septiembre-noviembre). Con esto confirmamos que la serie tiene un patrón de calendario: la gente viaja más en vacaciones y fiestas, menos en meses laborales.

2.2 Punto B

demanda_des_Air <- decompose(AirPassengers, type = "additive")
demanda_des_Air <- AirPassengers - demanda_des_Air$seasonal
demanda2 <- data.frame("demanda_ob" = as.numeric(AirPassengers)) %>% 
 mutate(demanda_descomp = as.numeric(demanda_des_Air)) %>%
  mutate(Fecha = fechas <- seq(
  from = as.Date("1949-01-01"),
  to   = as.Date("1960-12-01"),
  by   = "month"
))
ts_plot(demanda2, slider = T,
      title = "Demanda mensual de pasajeros",
        Ytitle= "Demanda de pasajeros",
        Xtitle= "Año-Mes")

Notamos en la grafica que la desestacionalización redujo bastante el efecto calendario. La serie original que vemos de color azul tiene picos fuertes cada año. Y en la serie desestacionalizada que observamos de color naranja suaviza esos picos. Y no se ven las subidas bruscas en verano ni las caídas en invierno. lo que queda es la tendencia de crecimiento más limpia. En resumidas cuentas deducimos que el efecto calendario se redujo mucho, casi desapareció, ahora la serie refleja sobre todo la tendencia de largo plazo y no los movimientos repetitivos de cada año.

2.3 Punto C

2.3.1 Correlograma con tres años de rezago

acf(AirPassengers, 
    main = "Correlograma ACF - Demanda Pasajeros", 
    lag.max = 36)

Notamos que el patron de correlacion disminuye de manera lenta, en el correlograma se ve que los valores de acf siguen altos incluso cuando aumenta el rezago. no caen de golpe, sino de manera gradual. eso pasa porque la serie tiene fuerte tendencia y estacionalidad, entonces cada observación está muy ligada a las anteriores durante mucho tiempo.

2.3.2 Prueba Dickey Fuller

dickey_full2 <- ur.df(AirPassengers, type = "trend", lags = 0)
summary(dickey_full2)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -74.244 -21.848   0.092  18.107 106.445 
## 
## Coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept) 25.95168    7.32595   3.542   0.000539 ***
## z.lag.1     -0.26819    0.05781  -4.639 0.00000795 ***
## tt           0.71076    0.16706   4.255 0.00003811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.65 on 140 degrees of freedom
## Multiple R-squared:  0.1333, Adjusted R-squared:  0.1209 
## F-statistic: 10.76 on 2 and 140 DF,  p-value: 0.00004483
## 
## 
## Value of test-statistic is: -4.6392 7.4143 10.764 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
dwtest(dickey_full2@testreg)
## 
##  Durbin-Watson test
## 
## data:  dickey_full2@testreg
## DW = 1.2541, p-value = 0.000001791
## alternative hypothesis: true autocorrelation is greater than 0
par(mfrow = c(1, 2))
acf(AirPassengers, )
pacf(AirPassengers)

serie_dif2 <- diff(AirPassengers, differences = 1)

si la serie no es estacionaria, primero transformo los datos con logaritmos para estabilizar la varianza, después aplico una diferenciación estacional (lag 12 si es mensual) para eliminar los patrones repetidos cada año y finalmente una diferenciación regular para quitar la tendencia, al terminar reviso con pruebas como adf y con los gráficos acf si la serie ya es estacionaria.

2.4 Punto D

dic_full_diff2 <- ur.df(serie_dif2, type = "trend", lags = 0)
options(scipen = 999)
summary(dic_full_diff2)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -106.951  -18.672   -4.909   22.711   67.829 
## 
## Coefficients:
##              Estimate Std. Error t value           Pr(>|t|)    
## (Intercept)  1.519755   5.489055   0.277              0.782    
## z.lag.1     -0.694109   0.081209  -8.547 0.0000000000000203 ***
## tt           0.001296   0.066542   0.019              0.984    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.5 on 139 degrees of freedom
## Multiple R-squared:  0.3446, Adjusted R-squared:  0.3352 
## F-statistic: 36.55 on 2 and 139 DF,  p-value: 0.0000000000001762
## 
## 
## Value of test-statistic is: -8.5472 24.3672 36.5465 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
par(mfcol = c(1, 2))
acf(serie_dif2)
pacf(serie_dif2)

los correlogramas de la serie estacionaria muestran un corte en la acf en el primer rezago y un pico en el rezago 12, mientras que la pacf no presenta cortes claros, esto sugiere un componente ma(1) y un ma(1) estacional con periodo 12, por lo que el modelo arima adecuado es un arima(0,1,1)

2.4.1 Estimar modelo arima (0,1,1)

  arima2 <- arima(AirPassengers, c(0, 1, 1), seasonal = list(order = c(0,1,1), perid = 12))
  stargazer(arima2, type = "HTML")

2.5 Punto E

demanda2 <- demanda2 %>% mutate(resid = arima2$residuals,
                              demanda_est = as.numeric(demanda_ob-resid)) %>% 
  select(demanda_ob, demanda_est, Fecha)
ts_plot(demanda2, slider = T,
      title = "Demanda mensual de pasajeros normal y ajustada",
        Ytitle= " Demanda de pasajeros",
        Xtitle= "Año-Mes")

La gráfica nos da a entender que el modelo ARIMA ajusta bien la demanda mensual de pasajeros, la serie estimada sigue de cerca la serie observada, replicando tanto la tendencia creciente como la estacionalidad anual con picos regulares. La diferencia entre ambas líneas es mínima, lo que nos muestra que los residuos son pequeños y no sistemáticos. por ende nos da a entender que el modelo captura correctamente el comportamiento de la serie original y que el ajuste es preciso.

par(mfcol = c(1, 2))
hist(arima2$residuals)
qqnorm(arima2$residuals)
qqline(arima2$residuals, col = "blue")

jarque.bera.test(arima2$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima2$residuals
## X-squared = 12.481, df = 2, p-value = 0.001949

El histograma está centrado en cero, pero el Noemal Q-Q Plot muestra que los residuos no siguen una distribución normal, ya que hay desviaciones importantes en las colas. esto indica que el modelo arima básico no cumple con el supuesto de normalidad.

ts_plot(demanda2, slider = T,
      title = "Demanda mensual de pasajeros normal y ajustada",
        Ytitle= "Demanda de pasajeros",
        Xtitle= "Año-Mes")

La gráfica muestra que el modelo ARIMA ajustado con auto.arima representa muy bien la demanda mensual de pasajeros, la demanda estimada sigue casi de forma exacta a la observada en toda la serie, incluyendo la tendencia creciente y los patrones estacionales. Las diferencias entre ambas curvas son mínimas, podemos indicar que los residuos son pequeños y que el modelo tiene un buen nivel de precisión para explicar y predecir la serie.

par(mfcol = c(1, 2))
hist(arima_mod2$residuals)
qqnorm(arima_mod2$residuals)
qqline(arima_mod2$residuals, col = "blue")

jarque.bera.test(arima_mod2$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima_mod2$residuals
## X-squared = 15.131, df = 2, p-value = 0.0005181
options(scipen = 999)

Muestra los residuos del modelo seleccionado por auto.arima. en este caso, los residuos siguen mejor la distribución normal, el histograma es más simétrico y en Normal Q-Q Plot los puntos están más alineados con la diagonal, aunque aún se observan ligeras desviaciones en los extremos.

2.5.1 Pronostico para los siguientes 5 meses

valores_fut2 <- forecast(arima_mod2, h = 5, level = c(0.95))
print(valores_fut2)
##          Point Forecast    Lo 95    Hi 95
## Jan 1961       445.6349 423.0851 468.1847
## Feb 1961       420.3950 393.9304 446.8596
## Mar 1961       449.1983 419.4892 478.9074
## Apr 1961       491.8399 460.0092 523.6707
## May 1961       503.3945 469.9953 536.7937
par(mfcol = c(1, 1))
plot(valores_fut2)

Muestra el pronóstico de la serie usando el modelo auto.arima. Se utilizó un arima(2,1,1)(0,1,0)[12], que incorpora diferenciación estacional. el pronóstico indica que el número de pasajeros seguirá aumentando, manteniendo el patrón de crecimiento y la estacionalidad observada en los años anteriores. los próximos 5 meses se proyectan con valores cercanos a 400–500 pasajeros.