Punto A

library(ggplot2)
library(readr)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
## 
## 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)
## Warning: package 'timetk' was built under R version 4.4.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.4.3
library(urca)
## Warning: package 'urca' was built under R version 4.4.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Cargando paquete requerido: zoo
## Warning: package 'zoo' was built under R version 4.4.3
## 
## 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)
## Warning: package 'tseries' was built under R version 4.4.3
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
energy <- read_excel("C:/Users/leona/OneDrive/Pictures/Desktop/Econometria 2/energy.xlsx")
getwd()
## [1] "C:/Users/leona/OneDrive/Pictures/Desktop/Econometria 2"
setwd("C:/Users/leona/OneDrive/Pictures/Desktop/Econometria 2")


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

ts_plot(serie_t, slider = T,
        color = "green",
        title = "Demanda mensual de energia en sincelejo 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")
ts_decompose(serie_t, type = "both")

Como podemos observar de la gráfica de la demanda real de energía en MWH de Sincelejo del mercado no regulado (enero-2000 hasta septiembre-2023). Así mismo, es difícil apreciar de cierta manera si existe o no patrones de estacionalidades claros, pero, se notan picos inusuales en ciertos lapsos de tiempo.

En el año 2002 hay un descenso rápido y es debido a un apagón debido a que la empresa que prestaba el servicio a Sincelejo, no solo a este sino a demás municipios y departamentos, Electro Caribe suspendió el pago de salarios atrasados a sus trabajadores, aunque dicen que no se debió a este suceso sino a un apagón que afectó a 18 municipios debido a la voladura de torres de energía.

Para el año 2020, se puede ver un aumento desproporcionado en la demanda de energía debido a la pandemia del covid-19, la cual paralizo al mundo entero haciendo que las personas tuvieran que estar confinadas en sus casas, lo que aumento el consumo residencial debido a que no se podía salir de sus casas hasta cierto tiempo, a su vez se liquidó a Electricaribe y el servicio fue asumido por Afinia y Air-e en la región caribe, en 2021 aumento debido a que se comenzó a medir y registrar mejor el servicio y además de la recuperación de la economía después de la pandemia.

#Punto B

demanda_descompensada <- decompose(serie_t, type = "additive")

demanda_des<- serie_t - demanda_descompensada$seasonal

demanda <- data.frame("demanda_ob" = as.numeric(serie_t)) %>% 
  mutate(demanda_des = as.numeric(demanda_des), Fecha = energy_m$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")

Como podemos observar al comparar ambas gráficas tanto la desestacionalizada y la normal, podemos analizar que el cambio fue mínimo ya que tienen mucha similitud en muchos puntos, por lo tanto, no se redujo tanto el efecto calendario.

#Punto C

df <- ur.df(serie_t, 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_t, )
pacf(serie_t)

serie_t_dif <- diff(serie_t, differences = 1)

La serie no es estacionaria al analizar a priori en el punto A, debido a que tiene una tendencia y no tiene unos patrones totalmente constantes sino que varían, también podemos observar en el correlograma en la AFC como va disminuyendo en el tiempo, pero muy lentamente y no es cercano a cero, por parte de la PAFC solo tiene un rezago significativo y los demás no lo son, es un AR(1) y el modelo ARIMA sería (1,0,0).

DF Según la prueba de hipótesis del nivel de significancia del coeficiente del test Dickey Fuller podemos decir que:

Ho = δ = 0 No es estacionaria Ho = δ ≠ 0 Es estacionaria

Podemos decir que según el resultado del coeficiente del rezago de la demanda mensual de energía (0.1693) este no es estadísticamente significativo, por lo cual la serie no es estacionaria. Por lo tanto, no se rechaza la hipótesis nula porque no existe suficiente evidencia estadística para afirmar lo contrario.

Según el test de DW nos indicia que no existe autocorrelación significativa entre los residuos, es decir, no hay evidencia estadística suficiente para rechazar la hipótesis nula de ausencia de autocorrelación.

Para volver una serie no estacionaria a una serie estacionaria solo se debe diferenciar una vez, en caso de que con una no basta se vuelve a diferenciar por segunda vez, hasta que se vuelva estacionaria.

#Punto D

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

  #Estimar modelo arima (0,1,1)
  
  arima_ma1 <- arima(serie_t, c(0, 1, 1), seasonal = list(order = c(0,1,1), perid = 12))
  
stargazer(arima_ma1, type = "html")
Dependent variable:
serie_t
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

Después de diferenciar la serie 1 sola vez, el modelo ARIMA (0,1,1), esta serie es estacionaria con solo una diferenciación (d = 1), podemos ver que en AFC nos muestra una media móvil de orden 1 (q = 1) mientras que PAFC no muestra ningún indicio de autoregresión (p = 0), este modelo es bueno, ya que no está sobre ajustado, lo cual nos permitirá tener un mejor pronostico y al mismo tiempo un buen análisis.

#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 energia en sincelejo normal y ajustada 2000-2023",
        Ytitle= "Demanda Real (MWh)",
        Xtitle= "Año")
par(mfcol = c(1, 2))
hist(arima_ma1$residuals)
qqnorm(arima_ma1$residuals)
qqline(arima_ma1$residuals, col = "blue")

jarque.bera.test(arima_ma1$residuals)
## 
##  Jarque Bera Test
## 
## data:  arima_ma1$residuals
## X-squared = 6658.9, df = 2, p-value < 0.00000000000000022
options(scipen = 99999)
arima_mod <- auto.arima(serie_t, 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")
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)
valores_fut <- forecast(arima_ma1, h = 5, level = c(0.95))
print(valores_fut)
##          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(valores_fut)

Modelo ARIMA (0,1,1)

Podemos observar en este modelo que gráficamente se nota que los errores o residuos siguen una distribución normal, pero en el test de Jarque-Bera el valor es de (6658.9) lo que nos dice que hay una desviación de la normalidad y su p-valor (0.00000000000000022) es menor a 0.05, por lo tanto, los residuos en el modelo ARIMA (0,1,1) no siguen una distribución normal, así mismo hay evidencia estadística suficiente para rechazar la hipótesis nula de que los residuos siguen una distribución normal.

Auto-ARIMA (0,1,0)

El auto arima determino que el mejor modelo no tiene ningún componente autorregresivo (p = 0), la serie esta diferencia 1 una vez (d = 1) y no tiene una media móvil (q = 0), ahora si nos vamos al gráfico podemos ver como la serie estimada tiene mucha similitud con la serie y gráfica original. Además, gráficamente se puede ver que los residuos tienen normalidad, pero el test de Jarque-Bera dio un resultado de (9670,1) un valor demasiado grande y su p-valor (0.00000000000000022) por lo tanto, los residuos en el modelo ARIMA (0,1,0) no siguen una distribución normal, así mismo hay evidencia estadística suficiente para rechazar la hipótesis nula de que los residuos siguen una distribución normal. No obstante, analizando un poco a detalles el histograma y el Q-Q Plot de los residuos podemos ver la existencia de valores atípicos tanto en el histograma como en el Q-Q Plot.

Pronostico para el modelo ARIMA (0,1,1)

Como se mostró en los análisis anteriores de los dos modelos, donde sus residuos no presentaron normalidad, se usará el modelo ARIMA (0,1,1), debido a que no muestra un sobre ajuste. Así mismo, se va a pronosticar los valores de los próximos 5 meses con este modelo, al realizar el pronóstico nos damos cuenta de que tiene valores estables entre 226 y 230 de MWh, sin mostrar si va a decrer o incrementar el consumo de esta. Sin embargo, el mes de octubre se puede esperar un consumo de 227.8084, con un lapso o espacio de 214.6709 MWh a 240.9459 MWh, según los valores dados, se espera que se mantenga en los próximos meses noviembre y diciembre, y que pueda llevar a subir el precio para enero y febrero respectivamente. Por lo tanto, el pronóstico sugiere que, aunque la demanda se mantiene relativamente estable, se espera un ligero repunte en 2024, lo cual puede generar presiones en los precios de la electricidad en Sincelejo si la oferta no se ajusta proporcionalmente.