Predicción Demanda Energía Eléctrica Hora 11
Aymara
13/4/2022
Presentación del caso
El objetivo principal de este trabajo es poder realizar predicciones a uno o dos días vista utilizando un modelo de predicción óptimo que nos asegure un margen de error lo más mínimo posible. Para ello utilizaremos un modelo de serires temporales teniendo en cuenta distintos factores como la temperatura o labolaridad y nos centraremos, en nuestro caso, en el tramo horario de las 11:00.
Lectura y preparación de datos
Antes del análisis vamos a llevar a cabo una lectura y preparación de los datos. Para ello, primero cargaremos en R las librerías con las que vamos a trabajar.
library(R.matlab)## Warning: package 'R.matlab' was built under R version 4.1.2## R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.##
## Attaching package: 'R.matlab'## The following objects are masked from 'package:base':
##
## getOption, isOpenlibrary(xts)## Warning: package 'xts' was built under R version 4.1.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.1.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(lubridate)## Warning: package 'lubridate' was built under R version 4.1.2##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(forecast)## Warning: package 'forecast' was built under R version 4.1.2## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(dygraphs)## Warning: package 'dygraphs' was built under R version 4.1.2Por otro lado, nuestros archivos con los datos son dos, uno con los datos de la demanda eléctrica para cada franja horaria y otro con las temperaturas.
dat0 = readMat("Demanda_PENR2.mat") #el archivo donde están los datos de demanda
dat = as.data.frame(dat0$Demanda)
colnames(dat)=c("fecha","hora","dem","nose")Para el procesado de los datos de demanda es necesario corregir el cambio horario que se da durante los solsticios de verano e invierno. Para los días que se retrasa la hora aparece que la demanda en esa hora es 0, por eso le asignamos el valor de la hora siguiente. En cambio, los días que se adelanta la hora hay el doble de demanda por eso la dividimos entre dos.
sel = which(dat$dem ==0)
dat$dem[sel] = dat$dem[sel+1]
i=which(diff(dat$dem)>10000)
dat$dem[i+1] = dat$dem[i+1]/2## Por horas
f0 = dat$fecha[dat$hora == 1] #poner en forma matricial las horas
n = length(f0)
dem = matrix(0,n-1,25)
dem[,1] = f0[-n]
for (h in 1:24){
x = dat$dem[dat$hora == h]
dem[(1:(n-1)),h+1] = x[1:(n-1)]
}# el ultimo dia no esta completo, cojo n-1 dias
dem = as.data.frame(dem)
names(dem) = c("fecha",paste0("h",1:24))Para realizar nuestra estimación del modelo nos encontramos con la problemática del covid19, ya que durante los meses de encierro la demanda energética fue menor porque salvo los servicios esenciales estaba todo cerrado. Por eso para nuestro estudio vamos a usar los datos de los años 2011 al 2018 para estimar nuestro modelo y el año 2019 para testearlo. No obstante, después de validar nuestro modelo en el año 2019 lo probaremos en el 2020 usando los datos del 2012 al 2019, y también lo probaremos en el 2021 usando los datos del 2013 al 2020.
## elegimos de 2011 a 2021
inicio = which(dem$fecha==20110101) #inicio estimacion
fin = which(dem$fecha==20211231) # fin prediccion
dem = dem[inicio:fin,] #nuestro archivo de demandaFactores que influyen de la demanda energética
Temperatura
La temperatura es un factor muy importante a tener en cuenta ya que parece razonable que en los casos extremos de mucho frío o mucho calor, la demanda crezca por el uso de aparatos de calefacción o refrigeración, respectivamente.
Por tanto, se distinguen tres categorías de meses, de noviembre a abril que muestran sensibilidad al frío, de junio a septiembre que muestran sensibilidad al calor y el resto serían los meses que corresponden con los climas templados del otoño y la primavera.
La sensibilidad a la temperatura no es constante a lo largo del tiempo, sino que se ve modificada por cambios en los hábitos de los consumidores, en el nivel de renta… que permitirá compra o no de equipamientos electrónicos y por innovaciones tecnológicas que obtengan una mejor eficiencia energética.
## Construyo la matriz de regresores
## desde 20110101
## hasta 20191231
## Temperaturas
dat0 = readMat("TempPrev_PenR2.mat")
fecha=dat0$TempPrev[[1]][,1]
tmax=rowMeans( dat0$TempPrev[[1]][,2:11] ) #calculo la media
tmin=rowMeans( dat0$TempPrev[[2]][,2:11] )
xmat = data.frame(fecha,tmax,tmin)
plot(xmat$tmax,type="l",col="darkseagreen2",ylab="Temperatura")
lines(xmat$tmin,col="hotpink3")
Labolaridad
Generalmente de lunes a viernes el cosumo es mayor, disminuyendo los fines de semana. Adicionalmente, la existencia de festivos nacionales o regionales afecta mucho a la evolución de la demanda, ya que el consumo eléctrico disminuye considerablemente.
En España el número de días festivos es considerable y representa aproximadamente el 10% de los días del año. Esto supone un aumento en la complejidad a la hora de modelar la demanda de energía eléctrica la cual ya de por sí resulta complicada para los días laborables.
El efecto de las festividades no sólo influye en el día en concreto en el que se produce si no que repercute en la demanda de energía eléctrica de los días anteriores y posteriores. Además, esta incidencia en los días contiguos el festivo varía en función del día de la semana en el que se produzca.
## selecciono los festivos entre
fiesta = readMat("Fiestas_PENR2.mat")
fiestas = fiesta$Fiestas #el peso de la población
i = which(fiestas[,1]>=20100924 & fiestas[,1]<=20220210)
fiestas = fiestas[i,]
nfest = dim(fiestas)[1]
xmat$fest = 0
for (k in 1:nfest){
j = which(fiestas[k,1] == xmat$fecha)
xmat$fest[j] = fiestas[k,2]
}
inicio2 = which(xmat$fecha==20110101) #inicio estimacion temp
fin2 = which(xmat$fecha==20211231) # fin prediccion temp
xmat = xmat[inicio2:fin2,]Demanda en horas previas
Vamos a trabajar con la hora 11, no obstante, al preparar los datos vamos a escoger también las 9:00 y 10:00 ya que más adelante veremos que la demanda energética en cierta hora está relacionada con las horas anteriores. Además vamos a dejar ya estandarizados los valores de la demanda en las horas previas.
datos = cbind(xmat,demand=dem[,10:12])
#la hora con la que trabajamos es la 11
#pero vamos a seleccionar también los datos de las 2 horas anteriores
View(datos) #hasta aquí hemos preparado nuestros datos desde 2011 hasta 2019
n = which(datos$fecha == 20181231) #fin de estimación
in12 = which(datos$fecha == 20120101) #inicio año 2012
in13 = which(datos$fecha == 20130101) #inicio año 2013
fin19 = which(datos$fecha == 20191231) #fin 2019
fin20 = which(datos$fecha == 20201231) #fin 2020
fin21 = which(datos$fecha == 20211231) #fin 2021
y = ts(datos,frequency = 7) #declaro como variables temporales
#estandarizo la demanda de las horas 9 y 10
h9max = max(datos$demand.h9); h9min = min(datos$demand.h9);
h10max = max(datos$demand.h10); h10min = min(datos$demand.h10);
z9 = (y[,"demand.h9"] - h9min)/(h9max-h9min)
z10 = (y[,"demand.h10"]-h10min)/(h10max-h10min)Modelos
Modelo 0 - Auto.Arima
Para estimar un primer Modelo 0 vamos a utilizar la función auto.arima incluida en la librería ‘forecast’ indicando que vamos a usar el logaritmo de la demanda para estabilizar la varianza.
m0 = auto.arima(y[,"demand.h11"],approximation = FALSE,stepwise = FALSE,lambda = 0)Con esto obtenemos una serie ARIMA con 4 parámetros MA, un sMA y una diferencia estacional de periodo 7.
m0## Series: y[, "demand.h11"]
## ARIMA(0,0,4)(0,1,1)[7]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ma1 ma2 ma3 ma4 sma1
## 0.6433 0.4083 0.3133 0.1179 -0.8038
## s.e. 0.0166 0.0195 0.0180 0.0147 0.0189
##
## sigma^2 estimated as 0.003218: log likelihood=5817.07
## AIC=-11622.15 AICc=-11622.12 BIC=-11584.36Para escoger un buen modelo nos interesa que los residuos sean ruido blanco. Para comprobarlo ponemos mirar la ACF de los residuos y ver si las barras están dentro del intervalo marcado o también el p-valor, el cual nos interesa que sea mayor de 0.05.
checkresiduals(m0)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,4)(0,1,1)[7]
## Q* = 647.86, df = 9, p-value < 2.2e-16
##
## Model df: 5. Total lags used: 14AICc0=m0$aicc
error0=sqrt(m0$sigma2)*100
error0## [1] 5.672557En este caso se ve claramente que el p-valor es demasiado pequeño y a parte existen muchos valores atípicos, ya que todavía no hemos tenido en cuenta la influencia de la temperatura y la laboralidad, las cuales serán muy importantes. A parte, el error sería de un 5.69%
Antes de tener en cuenta los regresores de temperatura y laboralidad, vamos a ajustar los parámetros que hemos obtenido con la función auto.arima. Para ello vamos a apoyarnos en la ACF y PACF.
tsdisplay(diff(log(y[1:n,"demand.h11"]),lag=7)) #de seguro una diferencia estacional
En la ACF vemos un decrecimiento exponencial con el que podríamos suponer que la serie componentes autorregresivos. Si nos fijamos en la PACF podría haber 3 ó 6 AR, ya que la tercera raya se sale de los límites marcados, la 4ª y la 5ª no, pero la 6ª se vuelve a salir. En el caso de la 7ª no lo tendríamos en cuenta porque sería la parte estacional que se repite cada 7 días.
Vamos a estimar un nuevo modelo 0b que sería una ARIMA (6,0,4)(0,1,1)[7]. Claramente el p-valor ha mejorado que es lo que nos interesa. No obstante, esto puede ser debido a que hemos introducido demasiados parámetros. Para comproar si podemos obivar alguno podemos fijarnos en el t-valor (la estimación dividida entre el std error) y si nos da mayor de 2 es que es significativo.
m0b = Arima(y[1:n,"demand.h11"],order=c(6,0,4),seasonal=list(order=c(0,1,1),period=7),lambda=0)m0b## Series: y[1:n, "demand.h11"]
## ARIMA(6,0,4)(0,1,1)[7]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ma1 ma2
## 0.2295 0.9565 -0.2975 -0.3134 0.2643 0.0400 0.3353 -0.7306
## s.e. 0.3060 0.2008 0.1695 0.1191 0.0492 0.0612 0.3053 0.1442
## ma3 ma4 sma1
## -0.0086 0.3059 -0.9853
## s.e. 0.1236 0.1067 0.0032
##
## sigma^2 estimated as 0.002803: log likelihood=4423.84
## AIC=-8823.67 AICc=-8823.56 BIC=-8751.94checkresiduals(m0b)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(6,0,4)(0,1,1)[7]
## Q* = 27.667, df = 3, p-value = 4.267e-06
##
## Model df: 11. Total lags used: 14AICc0b=m0b$aicc
error0b=sqrt(m0b$sigma2)*100
error0b## [1] 5.294207Para los componentes autorregresivos regulares podemos no tener en cuenta el AR6 (0.04/0.0612<2), el AR5 (0.2643/0.0492>2) sí que sería significativo, así que tendremos en cuenta todos los componentes autorregresivos regulares hasta el 5. En el caso de los componentes regulares de media móvil el MA4 (0.3059/0.1067>2) sería significativo, así que no quitaremos ninguno.
m0c = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),lambda=0)m0c## Series: y[1:n, "demand.h11"]
## ARIMA(5,0,4)(0,1,1)[7]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 0.3504 0.8956 -0.3125 -0.2877 0.2541 0.2158 -0.7397 -0.0078
## s.e. 0.1563 0.1165 0.1565 0.0929 0.0489 0.1565 0.1216 0.1179
## ma4 sma1
## 0.2808 -0.9852
## s.e. 0.0841 0.0033
##
## sigma^2 estimated as 0.002803: log likelihood=4423.52
## AIC=-8825.03 AICc=-8824.94 BIC=-8759.28checkresiduals(m0c)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(5,0,4)(0,1,1)[7]
## Q* = 27.656, df = 3, p-value = 4.289e-06
##
## Model df: 10. Total lags used: 13AICc0c=m0c$aicc
error0c=sqrt(m0c$sigma2)*100
error0c## [1] 5.293922Con esto llegamos a un nuevo modelo 0c que es similar al 0b en cuanto al p-valor y al AICc pero nos hemos quitado un parámetro. Ya con este modelo vamos a ir introduciendo los regresores de temperatura y laboralidad, no obstante, más adelante podemos modificar los parámetros AR y MA para precisar más nuestra estimación.
Modelo 1 con Temperaturas
Para este modelo vamos a tener en cuenta el regresor temperatura. Para ello, utilizaremos los valores de las medias de temperatura máxima y mínima que calculamos durante el procesado de datos.
Los cuales vamos a estandarizar e incluiremos en la matriz de regresores X1 de tal forma que en dicha matriz quedarán reflejadas las medias de la temperatura máxima, mínima y el cuadrado de ambas.
w = weekdays(ymd(xmat$fecha)) #lubridate
sab = ( w == "sábado")
dom = (w == "domingo")
plot(xmat$tmax,dem$h11,col="paleturquoise2",xlab='Temperatura Máxima media', ylab='Demanda eléctrica en la hora 11')
points(xmat$tmax[sab],dem$h11[sab],col="olivedrab1")
points(xmat$tmax[dom],dem$h11[dom],col="lightpink2")
El introducir el cuadrado de la temperatura tiene sentido ya que si nos fijamos en el gráfico de la demanda eléctrica en función de la temperatura máxima se puede ver que sigue una distribución parabólica (si hiciésemos lo mismo con la temperatura mínima, obtendríamos un gráfico similar).
m = lm(dem$h11 ~ xmat$tmax + I(xmat$tmax^2))
summary(m)##
## Call:
## lm(formula = dem$h11 ~ xmat$tmax + I(xmat$tmax^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -13525 -2302 1200 2569 6848
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50053.636 699.171 71.59 <2e-16 ***
## xmat$tmax -1650.413 64.564 -25.56 <2e-16 ***
## I(xmat$tmax^2) 33.127 1.399 23.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3587 on 4015 degrees of freedom
## Multiple R-squared: 0.1682, Adjusted R-squared: 0.1678
## F-statistic: 406.1 on 2 and 4015 DF, p-value: < 2.2e-16A parte si realizamos una regresión lineal podemos comprobar que los regresores de temperatura y temperatura al cuadrado son significativos.
## Pongo temperaturas entre 0 y 1 por cuestiones numéricas. Estandarizo.
Tmaxx = max(xmat$tmax); Tminx = min(xmat$tmax);
Tmaxm = max(xmat$tmin); Tminm = min(xmat$tmin);
zmax = (y[,"tmax"] - Tminx)/(Tmaxx-Tminx)
zmin = (y[,"tmin"]-Tminm)/(Tmaxm-Tminm)
X1 = cbind(zmax, zmax^2, zmin, zmin^2)m1 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X1[1:n,]) m1## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 ma4
## 0.3639 0.8589 -0.3709 -0.2871 0.2744 0.1426 -0.7727 0.0733 0.3349
## s.e. 0.1235 0.1182 0.1126 0.0930 0.0382 0.1247 0.1126 0.0950 0.0788
## sma1 zmax zmax^2 zmin zmin^2
## -0.9842 -0.3763 0.2738 -0.3866 0.3858
## s.e. 0.0033 0.0608 0.0567 0.0589 0.0552
##
## sigma^2 estimated as 0.002645: log likelihood=4510.25
## AIC=-8990.5 AICc=-8990.34 BIC=-8900.84checkresiduals(m1)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 39.245, df = 3, p-value = 1.54e-08
##
## Model df: 14. Total lags used: 17AICc1=m1$aicc
error1=sqrt(m1$sigma2)*100
error1## [1] 5.142521Al incluir la temperatura hemos reducido el error a 5.14%, no obstante, los residuos todavía no son ruido blanco. Por eso vamos a probar a introducir la influencia de la temperatura máxima y mínima del día anterior junto con sus cuadrados.
zmax_ayer =zmax
zmin_ayer =zmin
for(k in 2:(length(zmax)))
{zmax_ayer[k] =zmax[k-1]
zmin_ayer[k] =zmin[k-1]}
X1b =cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2)
m1b = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X1b[1:n,])m1b## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 ma4
## 0.3989 0.8406 -0.4212 -0.2930 0.2823 0.0921 -0.7796 0.1221 0.3545
## s.e. 0.1180 0.1242 0.1037 0.0976 0.0375 0.1196 0.1146 0.0913 0.0812
## sma1 zmax zmax^2 zmin zmin^2 zmax_ayer zmax_ayer^2
## -0.9836 -0.3088 0.2338 -0.1693 0.2101 -0.3146 0.2481
## s.e. 0.0035 0.0617 0.0572 0.0680 0.0689 0.0701 0.0675
## zmin_ayer zmin_ayer^2
## -0.1851 0.1440
## s.e. 0.0594 0.0563
##
## sigma^2 estimated as 0.002596: log likelihood=4539.3
## AIC=-9040.6 AICc=-9040.34 BIC=-8927.03checkresiduals(m1b)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 47.036, df = 3, p-value = 3.415e-10
##
## Model df: 18. Total lags used: 21AICc1b=m1b$aicc
error1b=sqrt(m1b$sigma2)*100
error1b## [1] 5.095199El modelo ha mejorado, el error se ha reducido al igual que el AICc (-9040.34) aun así todavía podemos introducir como regresores la temperatura de dos días antes.
zmax_anteayer =zmax_ayer
zmin_anteayer =zmin_ayer
for(k in 2:(length(zmax_ayer)))
{zmax_anteayer[k] =zmax_ayer[k-1]
zmin_anteayer[k] =zmin_ayer[k-1]}
X1c =cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2)
m1c = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X1c[1:n,])m1c## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 ma4
## 0.4370 0.8529 -0.4670 -0.3188 0.2855 0.0489 -0.8163 0.1540 0.3841
## s.e. 0.1173 0.1267 0.1006 0.0993 0.0374 0.1188 0.1161 0.0899 0.0822
## sma1 zmax zmax^2 zmin zmin^2 zmax_ayer zmax_ayer^2
## -0.9830 -0.3062 0.2378 -0.2104 0.2431 -0.2386 0.1985
## s.e. 0.0035 0.0610 0.0568 0.0680 0.0692 0.0724 0.0694
## zmin_ayer zmin_ayer^2 zmax_anteayer zmax_anteayer^2 zmin_anteayer
## -0.0323 0.0526 -0.1929 0.1318 -0.1702
## s.e. 0.0692 0.0704 0.0705 0.0680 0.0588
## zmin_anteayer^2
## 0.1087
## s.e. 0.0562
##
## sigma^2 estimated as 0.002563: log likelihood=4560.33
## AIC=-9074.65 AICc=-9074.27 BIC=-8937.17checkresiduals(m1c)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 49.414, df = 3, p-value = 1.065e-10
##
## Model df: 22. Total lags used: 25AICc1c=m1c$aicc
error1c=sqrt(m1c$sigma2)*100
error1c## [1] 5.062171Con la temperatura de dos días antes hemos pasado de un error de 5.09% si solo tenemos en cuenta la temperatura del día antes, a un 5.06%. Evidentemente se ha reducido pero esta mejora ha sido muy pequeña. Po eso vamos a concluir aquí el estudio de la influencia de la temperatura y vamos a enfocarnos en el factor de la laboralidad.
Modelo 2 con los festivos
En esta primera estimación tendremos en cuenta los festivos de Año Nuevo, Día de Reyes, Día del Trabajador, Día de la Asunción, Día del Pilar, Día de Todos los Santos, Día de la Constitución, Navidad (25 Dic) y Nochevieja.
No obstante, vamos a dejar ya preparadas unas matrices denominadas “navidad” que incluye todos los días de las vacaciones de Navidad y “asunción” que incluye los días próximos a esa fecha. Ya que durante la Navidad la gran mayoría de la gente coge vacaciones y la festividad de la Asunción coincide con la segunda quincena de Agosto, que es otra fecha clave en la que las familias se van de vacaciones y más adelante comprobaremos que al estimar la demanda energética en esos días va a resultar más difícil.
fech24 = ymd(y[,"fecha"])
finde = as.numeric(weekdays(fech24)=="sábado" | weekdays(fech24)=="domingo")
labor = 1 - finde #para solo tener en cuenta cuando no caen en fin de semana
I_01Ene = as.numeric( day(fech24)==1 & month(fech24)==1)*labor
I_06Ene = as.numeric(day(fech24)==6 & month(fech24)==1)*labor
I_01May = as.numeric(day(fech24)==1 & month(fech24)==5)*labor
#Asunción
I_15Ago = as.numeric(day(fech24)==15 & month(fech24)==8)*labor
I_13Ago = as.numeric(day(fech24)==13 & month(fech24)==8)*labor
I_14Ago = as.numeric(day(fech24)==14 & month(fech24)==8)*labor
I_16Ago = as.numeric(day(fech24)==16 & month(fech24)==8)*labor
I_17Ago = as.numeric(day(fech24)==17 & month(fech24)==8)*labor
I_12Oct = as.numeric(day(fech24)==12 & month(fech24)==10)*labor
I_01Nov = as.numeric(day(fech24)==1 & month(fech24)==11)*labor
I_06Dic = as.numeric(day(fech24)==6 & month(fech24)==12)*labor
#Navidad
I_22Dic = as.numeric(day(fech24)==22 & month(fech24)==12)*labor
I_23Dic = as.numeric(day(fech24)==23 & month(fech24)==12)*labor
I_24Dic = as.numeric(day(fech24)==24 & month(fech24)==12)*labor
I_25Dic = as.numeric(day(fech24)==25 & month(fech24)==12)*labor
I_26Dic = as.numeric(day(fech24)==26 & month(fech24)==12)*labor
I_27Dic = as.numeric(day(fech24)==27 & month(fech24)==12)*labor
I_28Dic = as.numeric(day(fech24)==28 & month(fech24)==12)*labor
I_29Dic = as.numeric(day(fech24)==29 & month(fech24)==12)*labor
I_30Dic = as.numeric(day(fech24)==30 & month(fech24)==12)*labor
I_31Dic = as.numeric(day(fech24)==31 & month(fech24)==12)*labor
fest = cbind(I_01Ene,I_06Ene,I_01May,I_15Ago,I_12Oct,I_01Nov,
I_06Dic, I_25Dic, I_31Dic)
#solo incluimos el festivo del 15 de Agosto de la Asunción y el 25 y 31 de Diciembre de Navidad
navidad=cbind(I_22Dic,I_23Dic,I_24Dic,I_26Dic,
I_27Dic,I_28Dic,I_29Dic,I_30Dic) #resto de días de Navidad
asuncion=cbind(I_13Ago,I_14Ago,I_16Ago,I_17Ago)
#resto de días del puente de la Asunción
X2 = cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,fest)
m2 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X2[1:n,])m2## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 0.8110 -0.6616 0.8689 -0.1568 -0.0079 -0.1157 0.5156 -0.3242
## s.e. 0.3299 0.3678 0.2388 0.3325 0.1585 0.3297 0.1671 0.1441
## ma4 sma1 zmax zmax^2 zmin zmin^2 zmax_ayer
## -0.1503 -0.9775 -0.3188 0.2397 -0.1573 0.2207 -0.2498
## s.e. 0.1953 0.0038 0.0387 0.0358 0.0427 0.0435 0.0443
## zmax_ayer^2 zmin_ayer zmin_ayer^2 zmax_anteayer zmax_anteayer^2
## 0.2083 -0.0959 0.0840 -0.1397 0.1012
## s.e. 0.0426 0.0426 0.0433 0.0442 0.0426
## zmin_anteayer zmin_anteayer^2 fest.I_01Ene fest.I_06Ene fest.I_01May
## -0.1093 0.0580 -0.4436 -0.3066 -0.2395
## s.e. 0.0372 0.0355 0.0126 0.0106 0.0105
## fest.I_15Ago fest.I_12Oct fest.I_01Nov fest.I_06Dic fest.I_25Dic
## -0.2091 -0.2109 -0.2108 -0.1728 -0.3114
## s.e. 0.0096 0.0104 0.0104 0.0105 0.0108
## fest.I_31Dic
## -0.1042
## s.e. 0.0131
##
## sigma^2 estimated as 0.001005: log likelihood=5930.72
## AIC=-11797.44 AICc=-11796.71 BIC=-11606.16checkresiduals(m2)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 37.714, df = 3, p-value = 3.249e-08
##
## Model df: 31. Total lags used: 34AICc2=m2$aicc
error2=sqrt(m2$sigma2)*100
error2## [1] 3.169706Al incluir los festivos el error se ha reducido al 3.17% y el AICc (-11696.71) también ha disminuido. El p-valor ha aumentado, aunque todavía no se acerca al 0.05 para que podamos decir que los residuos son ruido blanco.
Modelo 3 con Semana Santa
Nuestro siguiente paso va a ser introducir en el modelo los días Jueves y Viernes Santo.
## semana santa
# 20110422
# 20120406
# 20130329
# 20140418
# 20150403
# 20160325
# 20170414
# 20180330
# 20190419
#los viernes santos
ss_v1 = as.numeric( day(fech24)==22 & month(fech24)==4 & year(fech24) == 2011)
ss_v2 = as.numeric( day(fech24)==6 & month(fech24)==4 & year(fech24) == 2012)
ss_v3 = as.numeric( day(fech24)==29 & month(fech24)==3 & year(fech24) == 2013)
ss_v4 = as.numeric( day(fech24)==18 & month(fech24)==4 & year(fech24) == 2014)
ss_v5 = as.numeric( day(fech24)==3 & month(fech24)==4 & year(fech24) == 2015)
ss_v6 = as.numeric( day(fech24)==25 & month(fech24)==3 & year(fech24) == 2016)
ss_v7 = as.numeric( day(fech24)==14 & month(fech24)==4 & year(fech24) == 2017)
ss_v8 = as.numeric( day(fech24)==30 & month(fech24)==3 & year(fech24) == 2018)
ss_v9 = as.numeric( day(fech24)==19 & month(fech24)==4 & year(fech24) == 2019)
ss_v = ss_v1+ ss_v2+ ss_v3 + ss_v4 + ss_v5 + ss_v6+ ss_v7+ ss_v8 + ss_v9
#los jueves santos
ss_j1 = as.numeric( day(fech24)==21 & month(fech24)==4 & year(fech24) == 2011)
ss_j2 = as.numeric( day(fech24)==5 & month(fech24)==4 & year(fech24) == 2012)
ss_j3 = as.numeric( day(fech24)==28 & month(fech24)==3 & year(fech24) == 2013)
ss_j4 = as.numeric( day(fech24)==17 & month(fech24)==4 & year(fech24) == 2014)
ss_j5 = as.numeric( day(fech24)==2 & month(fech24)==4 & year(fech24) == 2015)
ss_j6 = as.numeric( day(fech24)==24 & month(fech24)==3 & year(fech24) == 2016)
ss_j7 = as.numeric( day(fech24)==13 & month(fech24)==4 & year(fech24) == 2017)
ss_j8 = as.numeric( day(fech24)==29 & month(fech24)==3 & year(fech24) == 2018)
ss_j9 = as.numeric( day(fech24)==18 & month(fech24)==4 & year(fech24) == 2019)
ss_j = ss_j1+ ss_j2+ ss_j3 + ss_j4 + ss_j5 + ss_j6+ ss_j7+ ss_j8 + ss_j9
X3 = cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,fest,ss_v,ss_j)
m3 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X3[1:n,])m3## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 0.9070 -0.1254 -0.7280 1.0696 -0.2342 -0.3084 0.0788 0.7765
## s.e. 0.0527 0.0564 0.0233 0.0609 0.0550 0.0450 0.0341 0.0186
## ma4 sma1 zmax zmax^2 zmin zmin^2 zmax_ayer
## -0.6106 -0.9709 -0.3223 0.2510 -0.1713 0.2365 -0.2234
## s.e. 0.0535 0.0046 0.0340 0.0313 0.0377 0.0382 0.0390
## zmax_ayer^2 zmin_ayer zmin_ayer^2 zmax_anteayer zmax_anteayer^2
## 0.1870 -0.0935 0.0856 -0.1327 0.0922
## s.e. 0.0371 0.0375 0.0382 0.0388 0.0372
## zmin_anteayer zmin_anteayer^2 fest.I_01Ene fest.I_06Ene fest.I_01May
## -0.0899 0.0464 -0.4439 -0.3056 -0.2523
## s.e. 0.0327 0.0314 0.0113 0.0096 0.0099
## fest.I_15Ago fest.I_12Oct fest.I_01Nov fest.I_06Dic fest.I_25Dic
## -0.2126 -0.2193 -0.2192 -0.1658 -0.3245
## s.e. 0.0089 0.0096 0.0096 0.0098 0.0099
## fest.I_31Dic ss_v ss_j
## -0.0968 -0.2649 -0.1596
## s.e. 0.0119 0.0090 0.0089
##
## sigma^2 estimated as 0.0007874: log likelihood=6287.63
## AIC=-12507.27 AICc=-12506.44 BIC=-12304.03checkresiduals(m3)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 92.319, df = 3, p-value < 2.2e-16
##
## Model df: 33. Total lags used: 36AICc3=m3$aicc
error3=sqrt(m3$sigma2)*100
error3## [1] 2.805993Una vez más el modelo ha mejorado, ha disminuido el AICc (-12506.44) y el error ahora es de 2.81%. Sin embargo, todavía hay atípicos en los residuos.
Modelo 4 Más festivos
Vamos a introducir más festivos nacionales y locales.
xfest_n = as.numeric(datos$fest==1)
loc = (datos$fest>0) & (datos$fest<1)
xfest_l = 0*xfest_n
xfest_l[loc] = xmat$fest[loc]
X4 = cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,fest,ss_v,ss_j,xfest_n,xfest_l)
m4 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X4[1:n,])m4## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 0.3721 -0.1612 0.6786 0.1975 -0.1963 0.2808 0.4627 -0.2684
## s.e. 0.1744 0.1394 0.1130 0.1586 0.0565 0.1744 0.0915 0.1164
## ma4 sma1 zmax zmax^2 zmin zmin^2 zmax_ayer
## -0.2569 -0.9638 -0.3077 0.2323 -0.1410 0.1956 -0.2236
## s.e. 0.0835 0.0057 0.0279 0.0257 0.0307 0.0313 0.0318
## zmax_ayer^2 zmin_ayer zmin_ayer^2 zmax_anteayer zmax_anteayer^2
## 0.2023 -0.1059 0.0921 -0.1090 0.0824
## s.e. 0.0305 0.0305 0.0313 0.0318 0.0306
## zmin_anteayer zmin_anteayer^2 fest.I_01Ene fest.I_06Ene fest.I_01May
## -0.0818 0.0503 -0.3325 -0.2023 -0.1500
## s.e. 0.0268 0.0260 0.0103 0.0093 0.0092
## fest.I_15Ago fest.I_12Oct fest.I_01Nov fest.I_06Dic fest.I_25Dic
## -0.0997 -0.1054 -0.1031 -0.0749 -0.2465
## s.e. 0.0087 0.0092 0.0092 0.0087 0.0090
## fest.I_31Dic ss_v ss_j xfest_n xfest_l
## -0.1020 -0.1491 0.0234 -0.1131 -0.2284
## s.e. 0.0093 0.0089 0.0091 0.0051 0.0072
##
## sigma^2 estimated as 0.0005313: log likelihood=6863.42
## AIC=-13654.85 AICc=-13653.92 BIC=-13439.66checkresiduals(m4)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 105.99, df = 3, p-value < 2.2e-16
##
## Model df: 35. Total lags used: 38AICc4=m4$aicc
error4=sqrt(m4$sigma2)*100
error4## [1] 2.304994El modelo sigue mejorando. Tanto el AICc (-13653.92) como el error han descendido pero todavía queda algún atípico.
Modelo 5. Identificación de Atípicos
Vamos a identificar manualmente los atípicos que quedan.
out1 = as.numeric(abs(residuals(m4)) > .1)
out_i1=which(abs(residuals(m4)) > 0.1)
fech24[out_i1]## [1] "2011-01-08" "2011-08-15" "2011-10-31" "2011-12-25" "2012-01-01"
## [6] "2012-01-02" "2012-03-29" "2012-04-30" "2012-11-14" "2012-12-06"
## [11] "2012-12-24" "2012-12-25" "2013-01-06" "2014-05-02" "2014-12-08"
## [16] "2014-12-24" "2015-12-08" "2015-12-24" "2016-10-31" "2016-12-25"
## [21] "2017-01-01" "2017-01-02" "2018-04-30" "2018-12-08" "2018-12-24"Si nos fijamos, estos valores atípicos corresponden a días anteriores o posteriores a algún festivo, tiene sentido ya que seguramente la gente hizo puente durante esos días. Por tanto vamos a incluir estos atípicos en nuestra matriz de regresores y estimar de nuevo el modelo.
X5 = cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,fest,ss_v,ss_j,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,xfest_n,xfest_l,out1)
m5 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X5[1:n,])m5## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 1.0234 0.0159 -0.8997 0.8708 -0.0872 -0.3596 -0.0902 0.7987
## s.e. 0.1005 0.1456 0.1048 0.1153 0.0676 0.0976 0.1161 0.0716
## ma4 sma1 zmax zmax^2 zmin zmin^2 zmax_ayer
## -0.3522 -0.9552 -0.3021 0.2283 -0.1320 0.1904 -0.2288
## s.e. 0.0795 0.0072 0.0253 0.0234 0.0279 0.0285 0.0289
## zmax_ayer^2 zmin_ayer zmin_ayer^2 fest.I_01Ene fest.I_06Ene
## 0.2094 -0.1170 0.0984 -0.3469 -0.2140
## s.e. 0.0276 0.0276 0.0283 0.0093 0.0085
## fest.I_01May fest.I_15Ago fest.I_12Oct fest.I_01Nov fest.I_06Dic
## -0.1820 -0.1017 -0.1204 -0.1295 -0.0684
## s.e. 0.0085 0.0080 0.0084 0.0084 0.0079
## fest.I_25Dic fest.I_31Dic ss_v ss_j zmax_anteayer
## -0.2679 -0.1027 -0.1651 0.0289 -0.0752
## s.e. 0.0081 0.0085 0.0081 0.0083 0.0291
## zmax_anteayer^2 zmin_anteayer zmin_anteayer^2 xfest_n xfest_l
## 0.0532 -0.0993 0.0672 -0.0987 -0.2366
## s.e. 0.0279 0.0245 0.0237 0.0047 0.0065
## out1
## -0.0885
## s.e. 0.0037
##
## sigma^2 estimated as 0.0004455: log likelihood=7121.4
## AIC=-14168.79 AICc=-14167.82 BIC=-13947.62checkresiduals(m5)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Q* = 93.28, df = 3, p-value < 2.2e-16
##
## Model df: 36. Total lags used: 39AICc5=m5$aicc
error5=sqrt(m5$sigma2)*100
error5## [1] 2.110593Con este nuevo modelo podemos comprobar que el AICc ha seguido disminuyendo (-14167.82) al igual que el error que ahora es de 2.11%
Validación Modelo 5
Ahora vamos a validar nuestro modelo m5 con los datos reales que tenemos del 2019.
my = Arima(y[,"demand.h11"],xreg = X5,model = m5)
## cojo las predicciones del 2019
fcast = fitted(my)[(n+1):fin19] #estoy cogiendo las últimas
plot(y[(n+1):fin19,"demand.h11"],type="l",col="palevioletred2",lwd=2,xlab="Dia del año",ylab="Demanda eléctrica en MW") #datos reales
lines(fcast,col="gray9",lwd=1) #predicciones
title("Demanda eléctrica 2019")
legend("topright",legend = c("Demanda real","Predicción"),lwd =3,col=c("palevioletred2", "gray9"),bty ="n")
ti <- seq(from = as.POSIXct("2019-01-01"),
to = as.POSIXct("2019-12-31"), by = "day")
xt = data.frame(real=y[(n+1):fin19,"demand.h11"],pred=fcast)
xt = xts(xt,order.by = ti)
dygraph(xt) %>% dyRangeSelector()Casi todos los días tienen una buena estimación. Sin embargo, observamos que en las fechas próximas al 15 de Agosto (día de la Asunción) y en el mes de Diciembre, por Navidad, la predicción no tiene nada que ver con la realidad. Esto es debido a que durante esos periodos es cuando la gente se va más de vacaciones, por eso vamos a estimar un nuevo modelo m6 en el que incluiremos en nuestra matriz de regresores los días de Navidad y la Asunción.
Modelo 6. Navidad y Puente Asunción
X6 = cbind(zmax, zmax^2, zmin, zmin^2, zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,fest,ss_v,ss_j,xfest_n,xfest_l,out1,navidad,asuncion)
m6 = Arima(y[1:n,"demand.h11"],order=c(5,0,4),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X6[1:n,])m6## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(5,0,4)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): Se han producido NaNs## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 1.2857 -0.1074 -1.0640 0.9214 -0.0867 -0.6498 -0.1544 0.9034
## s.e. 0.1334 0.2277 0.2548 NaN NaN 0.1316 0.1339 0.1487
## ma4 sma1 zmax zmax^2 zmin zmin^2 zmax_ayer
## -0.3676 -0.9575 -0.3019 0.2284 -0.1317 0.1912 -0.2293
## s.e. NaN 0.0065 0.0246 0.0227 0.0272 0.0278 0.0283
## zmax_ayer^2 zmin_ayer zmin_ayer^2 zmax_anteayer zmax_anteayer^2
## 0.2069 -0.1218 0.1025 -0.0646 0.0451
## s.e. 0.0270 0.0270 0.0277 0.0281 0.0270
## zmin_anteayer zmin_anteayer^2 fest.I_01Ene fest.I_06Ene fest.I_01May
## -0.1022 0.0742 -0.3575 -0.2127 -0.1809
## s.e. 0.0237 0.0230 0.0092 0.0083 0.0084
## fest.I_15Ago fest.I_12Oct fest.I_01Nov fest.I_06Dic fest.I_25Dic
## -0.1275 -0.1196 -0.1286 -0.0703 -0.3067
## s.e. 0.0087 0.0083 0.0083 0.0077 0.0090
## fest.I_31Dic ss_v ss_j xfest_n xfest_l out1 navidad.I_22Dic
## -0.1206 -0.1677 0.0175 -0.1000 -0.2266 -0.0881 0.0113
## s.e. 0.0086 0.0077 0.0079 0.0046 0.0066 0.0039 0.0080
## navidad.I_23Dic navidad.I_24Dic navidad.I_26Dic navidad.I_27Dic
## 0.0004 -0.0442 -0.0719 -0.0361
## s.e. 0.0084 0.0091 0.0080 0.0083
## navidad.I_28Dic navidad.I_29Dic navidad.I_30Dic asuncion.I_13Ago
## -0.0286 -0.0294 -0.0376 -0.0219
## s.e. 0.0078 0.0085 0.0086 0.0082
## asuncion.I_14Ago asuncion.I_16Ago asuncion.I_17Ago
## -0.0421 -0.0391 -0.0272
## s.e. 0.0080 0.0080 0.0074
##
## sigma^2 estimated as 0.0004245: log likelihood=7197.93
## AIC=-14297.86 AICc=-14296.15 BIC=-14004.95El modelo ha mejorado pero ahora nos encontramos con una nueva problemática y es que nos aparecen NaN en algunos parámetros. Esto se debe a que la suma de los coeficientes AR es cercana a 1, lo que muestra que los parámetros están cerca del borde de la región de estacionariedad. Esto causará dificultades al tratar de calcular los errores estándar.
A parte llegamos a la conclusión de que hemos introducido coeficientes de más. Por eso quitaremos los AR4, AR5 y MA4 que son donde encontramos el problema.
m6 = Arima(y[1:n,"demand.h11"],order=c(3,0,3),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X6[1:n,])m6## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(3,0,3)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 sma1 zmax
## 0.1656 0.9095 -0.1446 0.4785 -0.4785 -0.0380 -0.9477 -0.3013
## s.e. 0.2067 0.0975 0.1955 0.2070 0.1474 0.0746 0.0096 0.0248
## zmax^2 zmin zmin^2 zmax_ayer zmax_ayer^2 zmin_ayer zmin_ayer^2
## 0.2277 -0.1334 0.1919 -0.2344 0.2120 -0.1239 0.1036
## s.e. 0.0229 0.0274 0.0280 0.0285 0.0272 0.0272 0.0279
## zmax_anteayer zmax_anteayer^2 zmin_anteayer zmin_anteayer^2
## -0.0673 0.0491 -0.0973 0.0707
## s.e. 0.0283 0.0272 0.0239 0.0232
## fest.I_01Ene fest.I_06Ene fest.I_01May fest.I_15Ago fest.I_12Oct
## -0.3582 -0.2141 -0.1826 -0.1272 -0.1199
## s.e. 0.0093 0.0084 0.0084 0.0087 0.0083
## fest.I_01Nov fest.I_06Dic fest.I_25Dic fest.I_31Dic ss_v ss_j
## -0.1292 -0.0709 -0.3071 -0.1224 -0.1672 0.0176
## s.e. 0.0083 0.0079 0.0091 0.0088 0.0080 0.0083
## xfest_n xfest_l out1 navidad.I_22Dic navidad.I_23Dic
## -0.0991 -0.2269 -0.0884 0.0103 -0.0028
## s.e. 0.0046 0.0067 0.0040 0.0080 0.0084
## navidad.I_24Dic navidad.I_26Dic navidad.I_27Dic navidad.I_28Dic
## -0.0441 -0.072 -0.0367 -0.0282
## s.e. 0.0092 0.008 0.0084 0.0082
## navidad.I_29Dic navidad.I_30Dic asuncion.I_13Ago asuncion.I_14Ago
## -0.0313 -0.0420 -0.0207 -0.0405
## s.e. 0.0088 0.0086 0.0083 0.0080
## asuncion.I_16Ago asuncion.I_17Ago
## -0.0380 -0.0263
## s.e. 0.0081 0.0075
##
## sigma^2 estimated as 0.0004268: log likelihood=7188.72
## AIC=-14285.45 AICc=-14283.94 BIC=-14010.48checkresiduals(m6)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(3,0,3)(0,1,1)[7] errors
## Q* = 137.54, df = 3, p-value < 2.2e-16
##
## Model df: 45. Total lags used: 48AICc6=m6$aicc
error6=sqrt(m6$sigma2)*100
error6## [1] 2.065954Tanto el error (2.07%) como el AICc (-14283.94) se han reducido. No obstante, los residuos aun no son ruido blanco.
Modelo 7. Influencia de la hora 9 y 10
Ya, para afinar nuestro modelo vamos a introducir en la matriz de regresores la influencia de las dos horas anteriores a nuestra hora de estudio.
X7 = cbind(zmax, zmax^2, zmin, zmin^2,
zmax_ayer, zmax_ayer^2, zmin_ayer, zmin_ayer^2,
zmax_anteayer, zmax_anteayer^2, zmin_anteayer, zmin_anteayer^2,
fest,ss_v,ss_j,xfest_n,xfest_l,out1,navidad,asuncion,z9,z10)
m7 = Arima(y[1:n,"demand.h11"],order=c(3,0,3),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X7[1:n,])m7## Series: y[1:n, "demand.h11"]
## Regression with ARIMA(3,0,3)(0,1,1)[7] errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 sma1 zmax
## -1.2369 -0.2010 0.3022 1.5643 0.6703 -0.1138 -0.4531 -0.0817
## s.e. 0.1139 0.1455 0.0643 0.1151 0.1744 0.0857 0.0204 0.0089
## zmax^2 zmin zmin^2 zmax_ayer zmax_ayer^2 zmin_ayer zmin_ayer^2
## 0.0667 0.0474 -0.0127 -0.0260 0.0230 -0.0172 0.0152
## s.e. 0.0083 0.0100 0.0102 0.0109 0.0104 0.0104 0.0105
## zmax_anteayer zmax_anteayer^2 zmin_anteayer zmin_anteayer^2
## 0.0167 -0.0199 -0.0245 0.0224
## s.e. 0.0102 0.0098 0.0085 0.0082
## fest.I_01Ene fest.I_06Ene fest.I_01May fest.I_15Ago fest.I_12Oct
## -0.0535 -0.0061 -0.0282 -0.0267 -0.0191
## s.e. 0.0043 0.0038 0.0035 0.0033 0.0034
## fest.I_01Nov fest.I_06Dic fest.I_25Dic fest.I_31Dic ss_v ss_j
## -0.0254 0.0239 -0.021 0.0042 -0.0287 -0.0106
## s.e. 0.0034 0.0032 0.004 0.0036 0.0032 0.0032
## xfest_n xfest_l out1 navidad.I_22Dic navidad.I_23Dic
## -0.0053 0.0019 -0.0106 -0.0064 -0.0033
## s.e. 0.0019 0.0031 0.0016 0.0032 0.0033
## navidad.I_24Dic navidad.I_26Dic navidad.I_27Dic navidad.I_28Dic
## 0.0063 0.0068 0.0065 0.0043
## s.e. 0.0036 0.0030 0.0031 0.0031
## navidad.I_29Dic navidad.I_30Dic asuncion.I_13Ago asuncion.I_14Ago
## 0.0006 0.0037 -0.0002 -0.0032
## s.e. 0.0033 0.0035 0.0031 0.0029
## asuncion.I_16Ago asuncion.I_17Ago z9 z10
## 0.0017 -0.0007 -0.3429 1.0239
## s.e. 0.0029 0.0028 0.0125 0.0135
##
## sigma^2 estimated as 6.547e-05: log likelihood=9933.25
## AIC=-19770.49 AICc=-19768.85 BIC=-19483.57checkresiduals(m7)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(3,0,3)(0,1,1)[7] errors
## Q* = 49.229, df = 3, p-value = 1.166e-10
##
## Model df: 47. Total lags used: 50AICc7=m7$aicc
error7=sqrt(m7$sigma2)*100
error7## [1] 0.8091219Hemos pasado de un error de 2.07% al 0.81%, lo cual es buenísimo. El AICc se ha reducido a -19768.85 y si nos fijamos en la ACF de los residuos, salvo una rayita son ruido blanco. Por tanto podemos concluir que este es el mejor modelo.
Coef_AICc=cbind(AICc0,AICc0b,AICc0c,AICc1,AICc1b,AICc1c,AICc2,AICc3,AICc4,AICc5,AICc6,AICc7)
errores=cbind(error0,error0b,error0c,error1,error1b,error1c, error2,error3,error4,error5,error6,error7)
Coef_AICc## AICc0 AICc0b AICc0c AICc1 AICc1b AICc1c AICc2
## [1,] -11622.12 -8823.563 -8824.944 -8990.339 -9040.339 -9074.273 -11796.71
## AICc3 AICc4 AICc5 AICc6 AICc7
## [1,] -12506.44 -13653.92 -14167.82 -14283.94 -19768.85errores## error0 error0b error0c error1 error1b error1c error2 error3
## [1,] 5.672557 5.294207 5.293922 5.142521 5.095199 5.062171 3.169706 2.805993
## error4 error5 error6 error7
## [1,] 2.304994 2.110593 2.065954 0.8091219Para comparar todos los modelos que hemos estimado tenemos todos los AICc y los errores. Se puede ver como van descendiendo. Excepto en el paso del Modelo 0 al Modelo 0b en el cual el AICc aumenta (el error si disminuye). Lo cual nos indica que seguramente habíamos sobredimensionado el modelo al incluir tantos coeficientes AR y AM, no obstante, en el modelo 6 nos dimos cuenta y eliminamos algunos parámetros.
Validación Modelo 7
Ahora vamos a validar este modelo 7, que claramente es el mejor, con los datos reales del 2019.
my = Arima(y[,"demand.h11"],xreg = X7,model = m7)
## cojo las predicciones del 2019
fcast = fitted(my)[(n+1):fin19] #estoy cogiendo las últimas
plot(y[(n+1):fin19,"demand.h11"],type="l",col="palevioletred2",lwd=2,xlab="Dia del año",ylab="Demanda eléctrica en MW") #datos reales
lines(fcast,col="gray9",lwd=1) #predicciones
title("Demanda eléctrica 2019")
legend("topright",legend = c("Demanda real","Predicción"),lwd =3,col=c("palevioletred2", "gray9"),bty ="n")
ti <- seq(from = as.POSIXct("2019-01-01"),
to = as.POSIXct("2019-12-31"), by = "day")
xt = data.frame(real=y[(n+1):fin19,"demand.h11"],pred=fcast)
xt = xts(xt,order.by = ti)
dygraph(xt) %>% dyRangeSelector()Efectivamente, ahora con este modelo corregido los valores reales coinciden con los predichos. Lo cual nos lleva a suponer que la adición de la demanda de las horas anteriores como regresores ha contribuido mucho a la mejora del modelo.
error = y[(n+1):fin19,"demand.h11"]-fcast
fech2019 = fech24[(n+1):fin19]
diasem = weekdays(fech2019)
diasem = factor(diasem,levels=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))
post = data.frame(fecha=fech2019,diasem = diasem, mes=month(fech2019), resi=error)
sem=tapply(post$resi, post$diasem, sd)
mes=tapply(post$resi, post$mes, sd)
barplot(mes,ylab ="Error en MW",col="plum1")
title("Error medio por mes 2019")
Ahora vamos analizar los errores medios por mes. Se puede observar que en el 2019 el mes más difícil de predecir son Enero y Noviembre. En cambio Septiembre es un mes que se ajusta muy bien al modelo.
barplot(sem,ylab ="Error en MW",col="seagreen2")
title("Error medio por dia de la semana 2019")
diames = tapply(post$resi, list(post$mes,post$diasem) , sd)Si analizamos el error cometido según el día de la semana, cuesta mucho predecir el lunes, seguramente por ser el día siguiente al fin de semana.
opar <- par(no.readonly = TRUE)
par(mar = c(5, 5, 4, 6))
matplot(diames,type="l",xlab ="mes",ylab ="Demanda eléctrica en MW",lty = c(1, 1, 1, 1, 1, 1, 1),lwd = c(2, 2, 2, 2, 2, 2, 2),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"))
title("Error medio por dia de la semana de cada mes 2019")
legend(x ="topright",inset = c(-0.27, 0), legend=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"), lty = c(1, 1, 1, 1, 1, 1, 1),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"),lwd =2,xpd = TRUE,bty ="n")
on.exit(par(opar))Al comparar los errores cometidos por mes y día. Los que más cuestan predecir son los lunes de Enero, los sábados de Abril y Noviembre y los jueves de Agosto. Por el contrario, es muy sencillo predecir los sábados de Septiembre y los viernes de Julio.
Predicción de la demanda en el 2020
Ahora, con los datos del 2012 al 2019 vamos a estimar la predicción de la demanda en el 2020 con nuestro modelo.
m2020 = Arima(y[in12:fin19,"demand.h11"],order=c(3,0,3),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X7[in12:fin19,])
my = Arima(y[,"demand.h11"],xreg = X7,model = m2020)
fcast = fitted(my)[(fin19+1):fin20] #estoy cogiendo las últimas
plot(y[(fin19+1):fin20,"demand.h11"],type="l",col="palevioletred2",lwd=2,xlab="Dia del año",ylab="Demanda eléctrica en MW") #datos reales
lines(fcast,col="gray9",lwd=1) #predicciones
title("Demanda eléctrica 2020")
legend("topright",legend = c("Demanda real","Predicción"),lwd =3,col=c("palevioletred2", "gray9"),bty ="n")
ti <- seq(from = as.POSIXct("2020-01-01"),
to = as.POSIXct("2020-12-31"), by = "day")
xt = data.frame(real=y[(fin19+1):fin20,"demand.h11"],pred=fcast)
xt = xts(xt,order.by = ti)
dygraph(xt) %>% dyRangeSelector()Este año es un poco complicado ya que durante mediados de Marzo y todo Abril fue cuando se produjo la cuarentena durante la pandemia del covid19. Esto se ve reflejado en nuestro modelo en el que la demanda desciende drásticamente durante esos meses ya que salvo los servicios indispensables estaba todo cerrado.
Podemos comprobar que nuestro modelo es muy bueno ya que se ajusto a los datos reales durante ese periodo.
error = y[(fin19+1):fin20,"demand.h11"]-fcast
fech2020 = fech24[(fin19+1):fin20]
diasem = weekdays(fech2020)
diasem = factor(diasem,levels=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))
post = data.frame(fecha=fech2020,diasem = diasem, mes=month(fech2020), resi=error)
sem=tapply(post$resi, post$diasem, sd)
mes=tapply(post$resi, post$mes, sd)
barplot(mes,ylab ="Error en MW",col="plum1")
title("Error medio por mes 2020")
Si nos fijamos en el error medio por mes. En el 2020 el mes de Abril, cuando se produjo el Estado de Alarma, es claramente el más difícil de predecir.
barplot(sem,ylab ="Error en MW",col="seagreen2")
title("Error medio por dia de la semana 2020")
diames = tapply(post$resi, list(post$mes,post$diasem) , sd)En la relación de error producido según el día. En este año es más complicado estimar los fines de semana respecto al resto de días.
opar <- par(no.readonly = TRUE)
par(mar = c(5, 5, 4, 6))
matplot(diames,type="l",xlab ="mes",ylab ="Demanda eléctrica en MW",lty = c(1, 1, 1, 1, 1, 1, 1),lwd = c(2, 2, 2, 2, 2, 2, 2),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"))
title("Error medio por dia de la semana de cada mes 2020")
legend(x ="topright",inset = c(-0.27, 0), legend=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"), lty = c(1, 1, 1, 1, 1, 1, 1),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"),lwd =2,xpd = TRUE,bty ="n")
on.exit(par(opar))Si comparamos los errores cometidos según mes y día los viernes y jueves de Abril y los domingos de Octubre son los más costosos de estimar frente a los viernes de Junio.
Predicción de la demanda en el 2020
m2021 = Arima(y[in13:fin20,"demand.h11"],order=c(3,0,3),seasonal=list(order=c(0,1,1),period=7),
lambda=0,xreg=X7[in13:fin20,])
my = Arima(y[,"demand.h11"],xreg = X7,model = m2021)
fcast = fitted(my)[(fin20+1):fin21] #estoy cogiendo las últimas
plot(y[(fin20+1):fin21,"demand.h11"],type="l",col="palevioletred2",lwd=2,xlab="Dia del año",ylab="MW") #datos reales
lines(fcast,col="gray9",lwd=1) #predicciones
title("Demanda eléctrica 2021")
legend("topright",legend = c("Demanda real","Predicción"),lwd =3,col=c("palevioletred2", "gray9"),bty ="n")
ti <- seq(from = as.POSIXct("2021-01-01"),
to = as.POSIXct("2021-12-31"), by = "day")
xt = data.frame(real=y[(fin20+1):fin21,"demand.h11"],pred=fcast)
xt = xts(xt,order.by = ti)
dygraph(xt) %>% dyRangeSelector()En el caso del 2021 la gráfica de la demanda se asemeja más a la forma de la del 2019, ya que, a pesar de que aun había (y hay) covid, ya en ese año se empezó a vacunar a la gente y a relajar las restricciones y abrir comercios de ocio en los cuales se consume mucha energía.
Se puede ver que la predicción que hemos hecho con nuestro modelo utilizando los datos del 2013 al 2020 se ajusta muy bien con la demanda real.
error = y[(fin20+1):fin21,"demand.h11"]-fcast
fech2021 = fech24[(fin20+1):fin21]
diasem = weekdays(fech2021)
diasem = factor(diasem,levels=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))
post = data.frame(fecha=fech2019,diasem = diasem, mes=month(fech2021), resi=error)
sem=tapply(post$resi, post$diasem, sd)
mes=tapply(post$resi, post$mes, sd)
barplot(mes,ylab ="Error en MW",col="plum1")
title("Error medio por mes 2021")
A la hora de estudiar el error medio por mes. Los meses de Enero y Diciembre (Navidad) han sido los más difíciles de estimar. Seguramente la dificultad de Enero se debe a Filomena, una nevada que asoló toda España y posiblemente se utilizase la calefacción mucho más que durante un mes normal.
barplot(sem,ylab ="Error en MW",col="seagreen2")
title("Error medio por dia de la semana 2021")
diames = tapply(post$resi, list(post$mes,post$diasem) , sd)Según los días de la semana no hay mucha diferencia entre predecir un día u otro.
opar <- par(no.readonly = TRUE)
par(mar = c(5, 5, 4, 6))
matplot(diames,type="l",xlab ="mes",ylab ="Demanda eléctrica en MW",lty = c(1, 1, 1, 1, 1, 1, 1),lwd = c(2, 2, 2, 2, 2, 2, 2),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"))
title("Error medio por dia de la semana de cada mes 2021")
legend(x ="topright",inset = c(-0.27, 0), legend=c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"), lty = c(1, 1, 1, 1, 1, 1, 1),col = c("lightsalmon3","hotpink1","dodgerblue","goldenrod1","sienna2","orchid3","palegreen2"),lwd =2,xpd = TRUE,bty ="n")
on.exit(par(opar))Si analizamos el error medio por día y mes casi todos los días de la semana de Enero y los jueves, viernes y sábado de Diciembre son los más difíciles de estimar mientras que los sábado de verano (Junio a Septiembre) son los que mejor se ajustan al modelo.
Predicción de la demanda en el 2022
Para finalizar esta sería una predicción de la demanda en el 2022. No obstante, no disponemos de datos reales para contrastar ya que es un año que aun no ha finalizado y existen factores como una pandemia o una nevada que podrían afectar notablemente a los resultados.
futuro2022 <- forecast(h=7,m2021, xreg = X7)
autoplot(futuro2022, include=100)
Conclusiones
Como conclusión podemos afirmar que lo que más influye en la demanda eléctrica a cierta hora depende mucho de la demanda de las horas previas ya que cuando hemos introducido esos regresores en el modelo ha sido cuando el error se ha reducido considerablemente.
No obstante, no se deben olvidar otros factores como la temperatura o la laboralidad. A parte, también sería interesante estudiar en un futuro los factores socioeconómicos ya que la demanda se puede ver muy afectada según el estilo de vida de la población.