Datos de energía eléctrica y temperatura

Se trata de datos aportados por una alumna del Master (Nerea Ruíz) en una edición precedente. Proporcionan consumos horarios de energía eléctrica y temperatura en una zona no especificada del norte de España. Concebiblemente, ambas series tienen relación debido a la esperable demanda residencial de electricidad para calefacción y refrigeración.

Se hace alguna manipulación preliminar, mostrando diversos modos de leer los datos, transformarlos, y en particular manejar periodos de muestreo intra-diarios, que no pueden ser gestionados con fechas de tipo Date.

Lectura de los datos

Nos convendrá utilizar un formato de serie temporal tal como el provisto pro zoo o su extensión xts. El paquete lubridate contiene también funciones útiles en el tratamiento de datos fechados.

setwd("/home/etptupaf/ft/docencia/MasterModMat/MaterialClase/Ejemplos")
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(xts)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

Lo más fácil es leer el fichero provisto kwhC.csv y añadir las fechas a continuación. Con fines didácticos, se hace de varios modos alternativos a continuación.

Veamos el formato de entrada (Se invocan recursos del sistema operativo, Linux: lo que sigue no funcionará en Windows).

system("head kwhC.csv", intern=TRUE)
##  [1] "Fecha,Hora,kWh,Temp"       "1999/01/01,1,1775471,-100"
##  [3] "1999/01/01,2,1677699,-100" "1999/01/01,3,1505854,-100"
##  [5] "1999/01/01,4,1347838,-100" "1999/01/01,5,1220755,-100"
##  [7] "1999/01/01,6,1126550,-100" "1999/01/01,7,1046618,-100"
##  [9] "1999/01/01,8,924080,-100"  "1999/01/01,9,884754,-100"

Hay una cabecera, la segunda columna es la hora, y la cuarta parece codificar los valores inexistentes como -100, de lo que nos habremos de ocupar.

Método 1. Probar directamente con read.zoo. Si las fechas están en la primera columna y en formato no muy esotérico, es fácil que podamos leerlas directamente.

temp    <- read.zoo(file="kwhC.csv",header=TRUE,sep=",")
## Warning in zoo(rval3, ix): some methods for "zoo" objects do not work if the
## index entries in 'order.by' are not unique
head(temp)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##            Hora     kWh Temp
## 1999-01-01    1 1775471 -100
## 1999-01-01    2 1677699 -100
## 1999-01-01    3 1505854 -100
## 1999-01-01    4 1347838 -100
## 1999-01-01    5 1220755 -100
## 1999-01-01    6 1126550 -100
str(temp)
## 'zoo' series from 1999-01-01 to 2001-04-30
##   Data: int [1:20423, 1:3] 1 2 3 4 5 6 7 8 9 10 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:3] "Hora" "kWh" "Temp"
##   Index:  Date[1:20423], format: "1999-01-01" "1999-01-01" "1999-01-01" "1999-01-01" "1999-01-01" ...

Ha hecho lo que ha podido; las fechas las ha leido (bien) en formato Date. Nos avisa de que hay fechas repetidas. No tiene modo de saber que la segunda columna debe formar parte de la fecha. La hora “1” es presumiblemente la primera del día, la que va del las 00:00 a las 01:00; le asignaremos el tiempo 00:30. Cambiaremos el formato de la fecha, la uniremos con la hora, y crearemos un “timestamp” con el que fechar las series:

temp.1  <- format(index(temp),format="%Y-%m-%d")
temp.2  <- sprintf("%02i:30:00",coredata(temp[,"Hora"])-1)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
fechas  <- strptime(paste(temp.1,temp.2),format="%Y-%m-%d %H:%M:%S")
elec    <- zoo(temp[,2:3],order.by=fechas)
## Warning in zoo(temp[, 2:3], order.by = fechas): some methods for "zoo" objects
## do not work if the index entries in 'order.by' are not unique

## Warning in zoo(temp[, 2:3], order.by = fechas): some methods for "zoo" objects
## do not work if the index entries in 'order.by' are not unique
head(elec)
##                         kWh Temp
## 1999-01-01 00:30:00 1775471 -100
## 1999-01-01 01:30:00 1677699 -100
## 1999-01-01 02:30:00 1505854 -100
## 1999-01-01 03:30:00 1347838 -100
## 1999-01-01 04:30:00 1220755 -100
## 1999-01-01 05:30:00 1126550 -100
tail(elec)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##                         kWh Temp
## 2001-04-30 20:30:00 2287974   11
## 2001-04-30 21:30:00 2528566   11
## 2001-04-30 22:30:00 2447241   10
## 2001-04-30 23:30:00 2215467    9
## <NA>                1711597   13
## <NA>                2016979   15
str(elec)
## 'zoo' series from 1999-01-01 00:30:00 to NA
##   Data: int [1:20423, 1:2] 1775471 1677699 1505854 1347838 1220755 1126550 1046618 924080 884754 875651 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:20423] "1999-01-01" "1999-01-01" "1999-01-01" "1999-01-01" ...
##   ..$ : chr [1:2] "kWh" "Temp"
##   Index:  POSIXlt[1:20423], format: "1999-01-01 00:30:00" "1999-01-01 01:30:00" "1999-01-01 02:30:00" ...

Se ha vuelto a quejar de valores repetidos, lo que no tendría ya que suceder. Investigando vemos que hay un NA en el paso del 31-10-1999 al 1-11-1999. Probablemente tiene que ver con el cambio de hora. En los datos originales aparece una hora “25”. Investiguemos los NA en las fechas y veremos lo que podemos hacer con ellas.

which(is.na(fechas))
## [1]  7296 16032
temp[7294:7298,]
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##            Hora     kWh Temp
## 1999-10-31   23 1822444   15
## 1999-10-31   24 1730952   17
## 1999-10-31   25 1711597   13
## 1999-11-01    1 1526415   14
## 1999-11-01    2 1350875   13
fechas[7294:7298]
## [1] "1999-10-31 22:30:00 CET" "1999-10-31 23:30:00 CET"
## [3] NA                        "1999-11-01 00:30:00 CET"
## [5] "1999-11-01 01:30:00 CET"
temp[16030:16034]
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##            Hora     kWh Temp
## 2000-10-29   23 2271206   13
## 2000-10-29   24 2132854   12
## 2000-10-29   25 2016979   15
## 2000-10-30    1 1772289   11
## 2000-10-30    2 1570716   11
fechas[16030:16034]
## [1] "2000-10-29 22:30:00 CET" "2000-10-29 23:30:00 CET"
## [3] NA                        "2000-10-30 00:30:00 CET"
## [5] "2000-10-30 01:30:00 CET"
fechas[7296] <- strptime("1999-11-01 00:15:00",format="%Y-%m-%d %H:%M:%S")
fechas[16032] <- strptime("2000-10-30 00:15:00",format="%Y-%m-%d %H:%M:%S")
elec <- zoo(temp[,2:3],order.by=fechas)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique

Método 2. Empalmar los componentes de la fecha con un ISOdate()

temp     <- read.zoo(file="kwhC.csv",header=TRUE,sep=",")
## Warning in zoo(rval3, ix): some methods for "zoo" objects do not work if the
## index entries in 'order.by' are not unique
temp.1   <- index(temp)
fechas   <- ISOdate(year(temp.1),month(temp.1),day(temp.1),coredata(temp[,"Hora"])-1,30,0)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
elec     <- zoo(temp[,2:3],order.by=fechas)
## Warning in zoo(temp[, 2:3], order.by = fechas): some methods for "zoo" objects
## do not work if the index entries in 'order.by' are not unique
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
head(elec)
##                         kWh Temp
## 1999-01-01 00:30:00 1775471 -100
## 1999-01-01 01:30:00 1677699 -100
## 1999-01-01 02:30:00 1505854 -100
## 1999-01-01 03:30:00 1347838 -100
## 1999-01-01 04:30:00 1220755 -100
## 1999-01-01 05:30:00 1126550 -100
tail(elec)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##                         kWh Temp
## 2001-04-30 20:30:00 2287974   11
## 2001-04-30 21:30:00 2528566   11
## 2001-04-30 22:30:00 2447241   10
## 2001-04-30 23:30:00 2215467    9
## <NA>                1711597   13
## <NA>                2016979   15
str(elec)
## 'zoo' series from 1999-01-01 00:30:00 to NA
##   Data: int [1:20423, 1:2] 1775471 1677699 1505854 1347838 1220755 1126550 1046618 924080 884754 875651 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:20423] "1999-01-01" "1999-01-01" "1999-01-01" "1999-01-01" ...
##   ..$ : chr [1:2] "kWh" "Temp"
##   Index:  POSIXct[1:20423], format: "1999-01-01 00:30:00" "1999-01-01 01:30:00" "1999-01-01 02:30:00" ...
temp[7290:7299,]
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique
##            Hora     kWh Temp
## 1999-10-31   19 1392581   19
## 1999-10-31   20 1631020   17
## 1999-10-31   21 1798346   16
## 1999-10-31   22 1832797   15
## 1999-10-31   23 1822444   15
## 1999-10-31   24 1730952   17
## 1999-10-31   25 1711597   13
## 1999-11-01    1 1526415   14
## 1999-11-01    2 1350875   13
## 1999-11-01    3 1252638   13
fechas[7290:7299]
##  [1] "1999-10-31 18:30:00 GMT" "1999-10-31 19:30:00 GMT"
##  [3] "1999-10-31 20:30:00 GMT" "1999-10-31 21:30:00 GMT"
##  [5] "1999-10-31 22:30:00 GMT" "1999-10-31 23:30:00 GMT"
##  [7] NA                        "1999-11-01 00:30:00 GMT"
##  [9] "1999-11-01 01:30:00 GMT" "1999-11-01 02:30:00 GMT"
fechas[7296] <- strptime("1999-11-01 00:15:00",format="%Y-%m-%d %H:%M:%S")
fechas[16032] <- strptime("2000-10-30 00:15:00",format="%Y-%m-%d %H:%M:%S")
elec <- zoo(temp[,2:3],order.by=fechas)
## Warning in zoo(rval, index(x)[i]): some methods for "zoo" objects do not work if
## the index entries in 'order.by' are not unique

Método 3. Generar fechas y añadirlas. También reemplazaremos los -100 por NA.

temp   <- read.csv(file="kwhC.csv",header=TRUE)
head(temp)
##        Fecha Hora     kWh Temp
## 1 1999/01/01    1 1775471 -100
## 2 1999/01/01    2 1677699 -100
## 3 1999/01/01    3 1505854 -100
## 4 1999/01/01    4 1347838 -100
## 5 1999/01/01    5 1220755 -100
## 6 1999/01/01    6 1126550 -100
tail(temp)
##            Fecha Hora     kWh Temp
## 20418 2001/04/30   19 2241689   14
## 20419 2001/04/30   20 2238351   12
## 20420 2001/04/30   21 2287974   11
## 20421 2001/04/30   22 2528566   11
## 20422 2001/04/30   23 2447241   10
## 20423 2001/04/30   24 2215467    9
str(temp)
## 'data.frame':    20423 obs. of  4 variables:
##  $ Fecha: chr  "1999/01/01" "1999/01/01" "1999/01/01" "1999/01/01" ...
##  $ Hora : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ kWh  : int  1775471 1677699 1505854 1347838 1220755 1126550 1046618 924080 884754 875651 ...
##  $ Temp : int  -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 ...
temp[,4] <- ifelse(temp[,4]==-100,NA,temp[,4])
fechas   <- seq(ISOdate(1998,12,31,23,30), by = "hour", length.out = 20423)

Las fechas así generadas tienen clase ‘POSIX’.

class(fechas)
## [1] "POSIXct" "POSIXt"
fechas[7290:7299]
##  [1] "1999-10-31 16:30:00 GMT" "1999-10-31 17:30:00 GMT"
##  [3] "1999-10-31 18:30:00 GMT" "1999-10-31 19:30:00 GMT"
##  [5] "1999-10-31 20:30:00 GMT" "1999-10-31 21:30:00 GMT"
##  [7] "1999-10-31 22:30:00 GMT" "1999-10-31 23:30:00 GMT"
##  [9] "1999-11-01 00:30:00 GMT" "1999-11-01 01:30:00 GMT"

Podemos ahora fechar con ellas las series que tenemos:

elec.zoo <- zoo(temp[,3:4],order.by=fechas)
head(elec.zoo)
##                         kWh Temp
## 1998-12-31 23:30:00 1775471   NA
## 1999-01-01 00:30:00 1677699   NA
## 1999-01-01 01:30:00 1505854   NA
## 1999-01-01 02:30:00 1347838   NA
## 1999-01-01 03:30:00 1220755   NA
## 1999-01-01 04:30:00 1126550   NA
tail(elec.zoo)
##                         kWh Temp
## 2001-04-30 16:30:00 2241689   14
## 2001-04-30 17:30:00 2238351   12
## 2001-04-30 18:30:00 2287974   11
## 2001-04-30 19:30:00 2528566   11
## 2001-04-30 20:30:00 2447241   10
## 2001-04-30 21:30:00 2215467    9
str(elec.zoo)
## 'zoo' series from 1998-12-31 23:30:00 to 2001-04-30 21:30:00
##   Data: int [1:20423, 1:2] 1775471 1677699 1505854 1347838 1220755 1126550 1046618 924080 884754 875651 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:2] "kWh" "Temp"
##   Index:  POSIXct[1:20423], format: "1998-12-31 23:30:00" "1999-01-01 00:30:00" "1999-01-01 01:30:00" ...

Convertir a otros tipos es simple:

datos <- coredata(elec.zoo)
head(datos)
##          kWh Temp
## [1,] 1775471   NA
## [2,] 1677699   NA
## [3,] 1505854   NA
## [4,] 1347838   NA
## [5,] 1220755   NA
## [6,] 1126550   NA
class(datos)
## [1] "matrix" "array"
elec.xts <- xts(elec.zoo)

Como vimos en sesiones anteriores, el fechado facilita el alineado de series, su representación gráfica, agregación, etc. Representación gráfica con correcto rotulado de ejes Los plots de un objeto zoo y xts son como vimos diferentes.

plot(elec.zoo)

plot(elec.xts)

Podemos agregar con comodidad:

fechas <- index(elec.zoo)
elec.day.zoo <- aggregate(elec.zoo, by=as.POSIXct(format(fechas,format="%Y-%m-%d")), FUN=sum)
plot(elec.day.zoo)

elec.month.zoo <- aggregate(elec.zoo, by=as.yearmon, FUN=sum)
plot(elec.month.zoo)

Podemos tomar diferencias…

delec.day.zoo <- diff(elec.day.zoo,lag=1)
plot(delec.day.zoo)

…y reemplazar valores perdidos.

elec.day.zoo.miss <- elec.day.zoo
elec.day.zoo.miss[150:200,1] <- NA
plot(elec.day.zoo.miss)

elec.day.zoo.miss[,1] <- na.locf(elec.day.zoo.miss[,1])
plot(elec.day.zoo.miss)

(¡No se quiere decir que esto sea un modo recomendable de operar! Puede valer para unos pocos valores sueltos, pero si hay una racha de NA, mucho mejor emplear un suavizado de Kalman para interpolarlos.) Podemos calcular medias móviles, medianas móviles…

elec.day.zoo <- aggregate(elec.zoo, by=as.POSIXct(format(fechas,format="%Y-%m-%d")), FUN=sum)
elec.day.mm7.zoo <- rollapply(elec.day.zoo,width=15,align="center",FUN=mean)
plot(elec.day.mm7.zoo)

elec.day.med.zoo <- rollmedian(x=elec.day.zoo[,1],k=15)
plot(elec.day.zoo[,1])
lines(elec.day.mm7.zoo[,1],col="red")
lines(elec.day.med.zoo,col="green")

Veamos como crear, por ejemplo, un modelo con datos diarios que incorpore estacionalidad semanal mediante dummies:

require(dlm)
## Loading required package: dlm
crearMod <- function(x) {
  mod <- dlmModPoly(1,dV=exp(x[1]),dW=exp(x[2]))  +
         dlmModSeas(7,dV=0,dW=c(exp(x[3]),rep(0,5)))
  return(mod)
}
ajuste.elecday.lls <- dlmMLE(elec.day.zoo[,1],parm=rep(0,3), build=crearMod, 
                             method="Nelder-Mead")
ajuste.elecday.lls
## $par
## [1] 19.866667 35.633333  9.133333
## 
## $value
## [1] 15185.75
## 
## $counts
## function gradient 
##       38       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
mod.elecday.lls <- crearMod(ajuste.elecday.lls$par)
res.elecday.lls <- dlmFilter(elec.day.zoo[,1],mod.elecday.lls)
plot(elec.day.zoo[,1])
lines(res.elecday.lls$f,col="red")