library(readxl) 
library(ggplot2)
library(readr)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(timetk)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TSstudio)
library(urca)
library(lmtest)
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(tseries)
getwd()
## [1] "C:/Users/USUARIO/Desktop/Econometria 2"
setwd("C:/Users/USUARIO/Desktop/Econometria 2")
energy <- read_excel("energy.xlsx")

PUNTO A

energy_mensual <- energy %>% 
  summarise_by_time(.date_var = Fecha, .by = "month",
                                        demanda_real = mean(demanda_real, .na.rm = TRUE))
 
serie_tiempo <- ts (energy_mensual$demanda_real, start = c(2000, 1), end = c(2023, 9), 
              frequency = 12)

ts_plot(serie_tiempo, slider = TRUE,
        title = "Demanda mensual de energía en Sincelejo 2000-2023",
        Ytitle = "Demanda Real (MWh)",
        Xtitle = "Año")

#analisis: Haciendo el respectivo análisis, podemos observar cómo ha evolucionado la solicitud de energía cada mes en Sincelejo desde el año 2000 hasta el 2023, notamos que hay un aumento constante, pero también se aprecian ciertos ritmos marcados por las estaciones, que tienen que ver con el tiempo o clima y las tradiciones de la zona: normalmente, el gasto sube en los meses de más calor, como marzo, abril, agosto y septiembre, porque se usan más los aires acondicionados y los ventiladores debido al clima cálido, además en diciembre y enero por las fiestas de Navidad y el movimiento de las tiendas; sin embargo, baja un poco en los meses de lluvia, como mayo-junio y octubre-noviembre, cuando el calor afloja un poco y los negocios van más lentos, así que estos cambios a lo largo del año se deben a una mezcla de clima, sociedad y fechas señaladas.

ts_decompose(serie_tiempo, type = "both")

PUNTO B

demanda_descomp <- decompose(serie_tiempo, type = "additive")

demanda_des<- serie_tiempo - demanda_descomp$seasonal

demanda <- data.frame("demanda_ob" = as.numeric(serie_tiempo)) %>% 
  mutate(demanda_des = as.numeric(demanda_des), Fecha = energy_mensual$Fecha)

ts_plot(demanda, slider = T,
      title = "Demanda mensual de energia en sincelejo normal y desestacionalizada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")

#analisis: Al desestacionalizar la serie, el efecto calendario se redujo de forma significativa, ya que las variaciones que se repetían cada año en meses específicos (como los aumentos en temporadas calurosas o festivas y las caídas en meses lluviosos) desaparecieron casi por completo; de este modo, la nueva serie refleja con mayor claridad la evolución real de la demanda y los cambios propios de la tendencia, mientras que los patrones estacionales quedaron aislados en un componente aparte sin distorsionar el análisis.

PUNTO C

df <- ur.df(serie_tiempo, 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
dwtest(df@testreg)
## 
##  Durbin-Watson test
## 
## data:  df@testreg
## DW = 2.0209, p-value = 0.5228
## alternative hypothesis: true autocorrelation is greater than 0
par(mfrow = c(1, 2))
acf(serie_tiempo, )
pacf(serie_tiempo)

serie_tiempo_diferenciada <- diff(serie_tiempo, differences = 1)

#analisis: La serie de demanda mensual de energía no es estacionaria, ya que presenta una clara tendencia creciente y patrones estacionales que hacen que su media varíe en el tiempo; esto se confirma en el correlograma, donde la autocorrelación disminuye lentamente, y en la prueba de Dickey-Fuller, cuyo estadístico (-1.37) no supera los valores críticos, por lo que no se rechaza la hipótesis de raíz unitaria; en consecuencia, para lograr estacionariedad es necesario aplicar al menos una diferencia de primer orden y, de ser necesario, una diferenciación estacional que elimine los efectos de calendario.

PUNTO D

df_diff <- ur.df(serie_tiempo_diferenciada, 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(serie_tiempo_diferenciada)
pacf(serie_tiempo_diferenciada)

#Estimar modelo arima (0,1,1)
  
  arima_maI <- arima(serie_tiempo, c(0, 1, 1), seasonal = list(order = c(0,1,1), perid = 12))
  
  stargazer(arima_maI, type = "HTML")

#analisis: Los correlogramas y las pruebas nos confirman que, tras aplicar una primera diferencia la serie se vuelve estacionaria (ADF en la serie diferenciada: τ = -17.09, rechazo claro de raíz unitaria), y el patrón de autocorrelaciones apunta a una dinámica MA de primer orden: la ACF corta bruscamente en el lag-1 mientras que la PACF no muestra picos significativos persistentes, lo que sugiere p = 0, d = 1, q = 1 para la parte no estacional; además aparecen señales de dependencia anual (rezagos en 12, 24, …) que justifican incorporar una componente estacional, y el comportamiento estacional es consistente con una estructura MA(1) estacional; por tanto, el modelo que mejor describe la serie según correlogramas y pruebas es ARIMA(0,1,1)(0,1,1)[12], que captura tanto la corrección de choque de corto plazo como la regularidad anual de la demanda.

PUNTO E

demanda <- demanda %>% mutate(resid = arima_maI$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")
par(mfcol = c(1, 2))
hist(arima_maI$residuals)
qqnorm(arima_maI$residuals)
qqline(arima_maI$residuals, col = "orange")

jarque.bera.test(arima_maI$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima_maI$residuals
## X-squared = 6658.9, df = 2, p-value < 0.00000000000000022
modelo_arima <- auto.arima(serie_tiempo, 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 energia en sincelejo normal y ajustada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")
par(mfcol = c(1, 2))
hist(modelo_arima$residuals)
qqnorm(modelo_arima$residuals)
qqline(modelo_arima$residuals, col = "orange")

jarque.bera.test(modelo_arima$residuals)
## 
##  Jarque Bera Test
## 
## data:  modelo_arima$residuals
## X-squared = 9670.1, df = 2, p-value < 0.00000000000000022
pronostic <- forecast(arima_maI, h = 5, level = c(0.95))
print(pronostic)
##          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(pronostic)

#analisis: Con el modelo ARIMA(0,1,1)(0,1,1)[12] logramos observar adecuadamente los patrones de la demanda mensual de energía en Sincelejo, pues reproduce tanto la tendencia creciente como los patrones estacionales asociados a los ciclos anuales de consumo (altas temperaturas, vacaciones y festividades). La comparación entre la serie observada y la ajustada confirma un buen ajuste, ya que ambas curvas se superponen de forma cercana. Sin embargo, al evaluar los residuos mediante el test de Jarque-Bera, Encontramos que no siguen una distribución normal, lo que indica que persisten asimetrías y colas pesadas. Por esta razón, se recurrió al procedimiento automático de auto.arima, que ajusta de manera más flexible la especificación del modelo. Aun cuando los residuos tampoco alcanzan plena normalidad, el modelo optimizado mejora la capacidad predictiva y permite generar pronósticos confiables. Las proyecciones para los próximos cinco meses sugieren que la demanda continuará aumentando de manera sostenida, manteniendo los efectos de calendario típicos del consumo eléctrico en la región.