PUNTO A
energy_mensual <- energy %>%
summarise_by_time(.date_var = Fecha, .by = "month",
demanda_real = mean(demanda_real, na.rm = TRUE))
ts_energy <- ts (energy_mensual$demanda_real, start = c(2000, 1),
end = c(2023, 9), frequency = 12)
ts_plot(ts_energy, slider = T,
title = "Demanda mensual de energía en Sincelejo 2000-2023",
Ytitle = "Demanda real",
Xtitle = "Año")
ts_decompose(ts_energy, type = "both")
demanda_descom <- decompose(ts_energy, type = "additive")
Al analizar el ciclo de la gráfica de la demanda real mensual de
energía en mwh de la ciudad de Sincelejo del mercado no regulado
(enero-2000 hasta septiembre-2023), no se logra observar a simple vista
presencia de estacionalidad de la serie, debido que no hay evidencia de
un patrón recurrente de comportamientos en períodos específicos, pero sí
picos anormales en ciertos periodos, a excepción de una pequeña
tendencia ascendente en el consumo de energía durante los meses de
verano, enero, febrero y marzo.
En el 2002 se observa un descenso debido a un apagón causado por una
falla de una subestación de energía, que dejó sin suministro eléctrico a
18 departamentos. En 2006 al 2007 se presenta un pico negativo
influenciado por los efectos climáticos del fenómeno de El niño, según
un reporte presentado por el banco de la república, en ese año la
prolongada sequía aumentó el costo de los servicios y como consecuencia
se generó una disminución en el consumo energético. (Mejía A., 2007,
#)
Posteriormente, en 2020, se registró un incremento abrupto en la
demanda de energía. Éste aumento estuvo relacionado con las medidas de
confinamiento y salubridad pública adoptadas durante la pandemia de
COVID-19 (OPS, 2020). Tras la crisis sanitaria y la reactivación
económica, el consumo de energía no retornó a los niveles previos a la
pandemia, sino que se mantuvo elevado pues las personas siguieron su
nueva dinámica de consumo adoptada en la pandemia.
PUNTO B
demanda_des <- ts_energy - demanda_descom$seasonal
demanda <- data.frame("demanda_ob" = as.numeric(ts_energy))%>%
mutate(demanda_des = as.numeric(demanda_des),
Fecha = energy_mensual$Fecha)
ts_plot(demanda, slider = T,
title = "Demanda mensual de energía en Sincelejo \n normal y desestacionalizada \n 2000-2023",
Ytitle = "Demanda real",
Xtitle = "Año")
Al aplicar el proceso de desestacionalización a la serie temporal y
analizar la gráfica desestacionalizada con la normal, observamos que no
se obtuvo mucha diferencia con la original, ya que las gráficas son
similares en muchos puntos e indica que no se obtuvo una reducción del
efecto calendario.
PUNTO C
df <- ur.df(ts_energy, type = "trend", lags = 0)
summary(df)
##
## ###############################################
## # 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
par(mfcol = c(1, 2))
acf(ts_energy, )
pacf(ts_energy)
dwtest(df@testreg)
##
## Durbin-Watson test
##
## data: df@testreg
## DW = 2.0209, p-value = 0.5228
## alternative hypothesis: true autocorrelation is greater than 0
# Se aplica primera diferencia para volver la serie estacionaria
ts_energy_diff <- diff(ts_energy, differences = 1)
Tras una revisión previa de manera a priori en el punto A, se
determinó que la serie no es estacionaria. observamos en el ACF que
existe una disminución a medida que se aleja en el tiempo pero no es
cercano a cero, por lo que la serie no es estacionaria. En el PACF vemos
que hay un solo rezago significativo mientras que los demás no lo
son.
Dickey Fuller
Según este test, por medio de la prueba de hipótesis de
significancia individual del coeficiente, tenemos que:
Ho = ઠ = 0 ó p = 1 No es estacionaria.
Ho = ઠ ≠ 0 ó p ≠ 1 Es estacionaria
En este caso el resultado del coeficiente de demanda de energía es
de 0.1693, no es significativo y por lo tanto no es estacionaria.
Para volver una serie no estacionaria a estacionaria se hace una
diferenciación, y si sigue dando no estacionaria, se diferencia por
segunda vez hasta que sea estacionaria.
PUNTO D
df_diff <- ur.df(ts_energy_diff, type = "trend", lags = 0)
options(scipen = 999)
summary(df_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(ts_energy_diff, )
pacf(ts_energy_diff)
# Estimar el modelo ARIMA(0,1,1)
arima_ma1 <- arima(ts_energy, c(0, 1, 1),
seasonal = list(order = c(0,1,1), period = 12))
stargazer(arima_ma1, type = "html")
|
|
|
|
Dependent variable:
|
|
|
|
|
|
ts_energy
|
|
|
|
ma1
|
0.045
|
|
|
(0.069)
|
|
|
|
|
sma1
|
-1.000***
|
|
|
(0.077)
|
|
|
|
|
|
|
Observations
|
272
|
|
Log Likelihood
|
-916.639
|
|
sigma2
|
43.057
|
|
Akaike Inf. Crit.
|
1,839.279
|
|
|
|
Note:
|
p<0.1; p<0.05;
p<0.01
|
El modelo ARIMA (0,1,1), resulta ser el más eficiente ya que muestra
la dinámica de la serie de tiempo, esta serie al ya ser diferenciada se
volvió estacionaria, esta necesito una sola diferenciación (d=1), además
el ACF nos muestra un comportamiento de media móvil de orden 1 (q=1) por
otro lado el PACF no muestra la necesidad de un componente
autorregresivo (p=0), así este modelo es el más parsimonioso y además,
evita el sobre ajuste, lo que nos permite una interpretación mejor y un
pronóstico más exacto de esta serie.
PUNTO E
demanda <- demanda %>%
mutate(resid = arima_ma1$residuals,
demanda_est = as.numeric(demanda_ob-resid))%>%
select(demanda_ob, demanda_est, Fecha)
ts_plot(demanda, slider = T,
title = "Demanda mensual de energía en Sincelejo \n normal y ajustada \n 2000-2023",
Ytitle = "Demanda real",
Xtitle = "Año")
par(mfcol = c(1, 2))
hist(arima_ma1$residuals)
qqnorm(arima_ma1$residuals)
qqline(arima_ma1$residuals, col = "red")
jarque.bera.test(arima_ma1$residuals)
##
## Jarque Bera Test
##
## data: arima_ma1$residuals
## X-squared = 6658.9, df = 2, p-value < 0.00000000000000022
# Version automatica
modelo_arima <- auto.arima(ts_energy, 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 = modelo_arima$residuals,
demanda_est = as.numeric(demanda_ob-resid))%>%
select(demanda_ob, demanda_est, Fecha)
ts_plot(demanda, slider = T,
title = "Demanda mensual de energía en Sincelejo \n normal y ajustada \n 2000-2023",
Ytitle = "Demanda real",
Xtitle = "Año")
par(mfcol = c(1, 2))
hist(modelo_arima$residuals)
qqnorm(modelo_arima$residuals)
qqline(modelo_arima$residuals, col = "red")
jarque.bera.test(modelo_arima$residuals)
##
## Jarque Bera Test
##
## data: modelo_arima$residuals
## X-squared = 9670.1, df = 2, p-value < 0.00000000000000022
# Pronosticamos con el modelo ARIMA(0,1,1)
futureval <- forecast(arima_ma1, h = 5, level = c(0.95))
print(futureval)
## Point Forecast Lo 95 Hi 95
## Oct 2023 227.8084 214.6709 240.9459
## Nov 2023 226.3940 207.3944 245.3936
## Dec 2023 227.4301 203.9913 250.8688
## Jan 2024 223.1286 195.9668 250.2905
## Feb 2024 230.3347 199.9072 260.7621
par(mfcol = c(1, 1))
plot(futureval)
Modelo ARIMA (0,1,1)
Estudiando la normalidad de los errores en este modelo obtuvimos
como resultado el valor estadístico Jarque-Bera es de 6658.9, lo que
indica una desviación significativa de la normalidad, pues este valor
está bastante alejado de cero.Además, su p-value (0.00000000000000022)
es inferior a 0.05, lo que indica que los residuos no siguen una
distribución normal.
El histograma muestra una alta concentración en el centro con
valores extremos, mientras que el gráfico Q-Q plot revela desviaciones
en las colas. Estos hallazgos confirman el test de Jarque-Bera,
rechazando la normalidad.
AUTO-ARIMA (0,1,0)
Al momento de aplicar el auto-arima, el modelo se ajusta un poco más
al gráfico, se nota que la serie estimada se asemeja más a la serie
original. En cuanto al test de Jarque Bera, los resultados del X-squared
(9670.1) un valor grande alejado de cero y su p-value
(0.00000000000000022) menor a 0,05, esto indica los errores no siguen
una distribución normal.
En el histograma se observa una gran cantidad residuos que se
encuentran muy cerca de cero, esto indica la existencia de valores
atípicos. En el Q-Q Plot los datos presentan distorsión en los extremos,
lo que indica que los residuos no se distribuyen de manera normal.
Pronóstico del modelo ARIMA (0,1,1):
Dado que ambos modelos no presentaron normalidad en los errores,
procedimos a realizar el pronóstico utilizando el modelo con menor
cantidad de rezagos, que fue el modelo ARIMA (0,1,1), éste modelo se usó
para pronosticar los valores de los próximos 5 meses. Éste muestra
valores estables que oscilan entre 226 y 230 de MWh, sin una tendencia
clara hacia el crecimiento o decrecimiento del consumo. En octubre de
2023, se espera un consumo de 227.8 MWh, con un intervalo entre 214.67
MWh y 240.95MWh, según el pronóstico se espera que este patrón se
mantenga en los meses de fin de año noviembre y diciembre, como también
se espera que la demanda de consumo tenga cierta tendencia al aumento a
principios del año 2024, por el inicio de la temporada de verano.