El objeto de series temporales

Las fechas se interpretan desde a EPOCH:

#EPOCH : 1 de Enero de 1970

Formas de obtener la fecha del día en que estás:

#HOY
Sys.Date()
## [1] "2018-11-21"

Debemos indicar siempre el formato de la fecha:

#año dos dígitos
as.Date("1/1/80", format="%m/%d/%y")
## [1] "1980-01-01"
#año cuatro dígitos
as.Date("1/1/1980", format="%m/%d/%Y")
## [1] "1980-01-01"

Si omitimos el parámetro format, siempre debe seguir el siguiente formato:

#yyyy-mm-dd o yyyy/mm/dd
as.Date("2018-01-06")
## [1] "2018-01-06"

Si los intorucieramos mal, se guardaría mal internamente:

as.Date("88/05/19")
## [1] "0088-05-19"

Si queremos transformar en horas la fecha que guardamos desde el momento de la epoch, usamos la función as.numeric:

##fecha -> numero
as.numeric(as.Date("1988/05/19"))
## [1] 6713
# Quantos días  han pasado hoy desde el EPOCH?
as.numeric(Sys.Date())
## [1] 17856

Podemos usar otras formas de separar la información, pero siempre debemos avisar el formato. La fecha se debe introudcir en el idioma que está en el sistema o cambiarlo con el parámetro locale:

#Nombre de los meses
as.Date("Ene 6, 2018", format="%b %d, %Y")
## [1] NA
as.Date("Enero 6, 18", format="%B %d, %y")
## [1] "2018-01-06"

Miremos como sería contar o restar días desde la EPOCH:

#Fechas desde días de EPOCH
dt <- 2018
class(dt) <- "Date" # indico que es una fecha
dt
## [1] "1975-07-12"
dt <- -2018
class(dt) <- "Date"
dt
## [1] "1964-06-23"

Se pueden crear fechas desde días de un punto dado:

#Fechas desde días de un punto dado
as.Date(2018, origin = as.Date("1988-05-19"))
## [1] "1993-11-27"
as.Date(-2018, origin = as.Date("1988-05-19"))
## [1] "1982-11-09"

Trabajamos con componentes de las fechas:

#Componentes de las fechas
dt <- as.Date(2018, origin = as.Date("1988-05-19"))
dt
## [1] "1993-11-27"
#Año en 4 dígitos
format(dt, "%Y")
## [1] "1993"
#Año como número en lugar de String
as.numeric(format(dt, "%Y"))
## [1] 1993
#Año en 2 dígitos
format(dt, "%y")
## [1] "93"
#Año como número en lugar de String
as.numeric(format(dt, "%y"))
## [1] 93
#Mes como String
format(dt, "%b") # mes  abreviado
## [1] "nov."
format(dt, "%B") # mes ampliado
## [1] "noviembre"
months(dt)
## [1] "noviembre"
weekdays(dt)
## [1] "sábado"
quarters(dt)
## [1] "Q4"
julian(dt) # mes con el calendario juliano (desde la EPOCH)
## [1] 8731
## attr(,"origin")
## [1] "1970-01-01"
julian(dt, origin = as.Date("1980/06/25")) # podemos establecer la fecha que no sea la EPOCH con el parametro "origin"
## [1] 4903
## attr(,"origin")
## [1] "1980-06-25"

Otras sintáxis importantes para tratar con fechas:

Operaciones y secuencias de fechas

Trabajaremos con la fecha “1/1/2001”

dt <- as.Date("1/1/2001", format = "%d/%m/%Y")
dt+100 # sumamos 100 días
## [1] "2001-04-11"
dt-100 # restamos 10 días
## [1] "2000-09-23"
dt+31 # suma 31 días
## [1] "2001-02-01"

Trabajando con dos fechas:

dt
## [1] "2001-01-01"
dt2 <- as.Date("2001/01/02")
dt2-dt # restamos las dos fechas (1 día)
## Time difference of 1 days
dt-dt2
## Time difference of -1 days
as.numeric(dt-dt) # podemos sacar el número exacto
## [1] 0
# que fecha vino antes?
dt<dt2
## [1] TRUE
dt==dt2
## [1] FALSE
dt2<dt
## [1] FALSE
seq(dt, dt+365, "month") # cuenta de mes en mes 
##  [1] "2001-01-01" "2001-02-01" "2001-03-01" "2001-04-01" "2001-05-01"
##  [6] "2001-06-01" "2001-07-01" "2001-08-01" "2001-09-01" "2001-10-01"
## [11] "2001-11-01" "2001-12-01" "2002-01-01"
seq(dt, as.Date("2001/01/10"), "day") # cuenta de día en día
##  [1] "2001-01-01" "2001-01-02" "2001-01-03" "2001-01-04" "2001-01-05"
##  [6] "2001-01-06" "2001-01-07" "2001-01-08" "2001-01-09" "2001-01-10"
seq(dt, dt+365, "2 months") # cuenta bimensual
## [1] "2001-01-01" "2001-03-01" "2001-05-01" "2001-07-01" "2001-09-01"
## [6] "2001-11-01" "2002-01-01"
# una secuencía a partir de una fecha dada
seq(from = dt, by = "4 months", length.out = 6)[3] 
## [1] "2001-09-01"
seq(from = dt, by = "3 weeks", length.out = 6)[4]
## [1] "2001-03-05"

Análisis preliminar de una serie temporal

Analizaremos los datos de Wallmart:

wmt <- read.csv("../DataSets/WMT.csv", stringsAsFactors = F)
plot(wmt$Adj.Close, type = "l")

Podemos ver el valor de cierre de las acciones de Wallmart. Vamos a ver las diferencias de un día al siguiente:

d <- diff(wmt$Adj.Close, lag = 2) # usamos este parametro para encontrar diferencias dentro de un mismo periodo (comparamos  una observación con la anterior)
head(d)
## [1]  0.112281  0.224563  0.044912 -0.247019 -0.606312  0.067372

Vamos a representar las diferencias:

plot(d, type="l")

Podemos ver la representación de los movimientos de un día a otro. REalizemos un histograma de los cambios diarios que ha habido, y buscar si se aproxima a una curva normal con el objetivo de buscar si los cambios tienen cierta tendencia.

hist(d, prob=T, 
     ylim = c(0,0.8), # ajsutamos la altura máxima del histograma
     breaks = 30, main = "Walmart Stocks", col = "green")# esclamos a probabilidades para ver una buena forma  del hisograma
lines(density(d), lwd=2)

Hay una distribución bastante normal, y no vemos ningún fuera de lo común.

Analizemos los datos a nivel mensual, con el que crearemos una serie temporal:

wmt.m <- read.csv("../DataSets/WMT-monthly.csv", stringsAsFactors = F)
wmt.m <- wmt.m[2:nrow(wmt.m),] # eliminamos la primera fila que es null
wmt.m$Date <- as.Date(wmt.m$Date) # combertimos a date
wmt.m.ts <- ts(wmt.m$Adj.Close)

d <- diff(as.numeric(wmt.m.ts)) # vamos a crear las diferencias respecto la serie  temporal creada
wmt.m.return <- d / lag(as.numeric(wmt.m.ts), k = -1) # k = -1 excuyenndo un día que queda fuera por la diferencia, la función lag nos sirve para calcuar diferencias en intervalos mayores. Si queremos calcuar diferencias de dos días en dos días debemos poner  lag = 2. Hemos analizado el retorno por stock, de un mes al siguiente
## Warning in d/lag(as.numeric(wmt.m.ts), k = -1): longitud de objeto mayor no
## es múltiplo de la longitud de uno menor
hist(wmt.m.return, prob = T, col = "blue")

El objeto dde serio temporal de R

En este caso cubriremos las series temporales con un solo dato como con múltiples datos. Cargamos los datos:

s <- read.csv("../DataSets/ts-example.csv")
head(s)
s.ts <- ts(s) # pasamos el array de números al objeto temporal:
class(s.ts)
## [1] "ts"
head(s.ts)
## [1]  51  56  37 101  66  63
plot(s.ts)

Creamos una serie temporal incluyendo el punto de origen

s.ts.a <- ts(s, start = 2001) 
head(s.ts.a)
## [1]  51  56  37 101  66  63
plot(s.ts.a)

Podemos añadir datos mensuales añadiendo la freqüencia:

s.ts.m <- ts(s, start = c(2001,1), frequency = 12) # indicamos cuando empieza y la frecuencia mesnual
s.ts.m
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2001  51  56  37 101  66  63  45  68  70 107  86 102
## 2002  90 102  79  95  95 101 128 109 139 119 124 116
## 2003 106 100 114 133 119 114 125 167 149 165 135 152
## 2004 155 167 169 192 170 180 175 207 164 204 180 203
## 2005 215 222 205 202 203 209 200 199 218 221 225 212
## 2006 250 219 242 241 267 249 253 242 251 279 298 260
## 2007 269 257 279 273 275 314 288 286 290 288 304 291
## 2008 314 290 312 319 334 307 315 321 339 348 323 342
## 2009 340 348 354 291
plot(s.ts.m)

Indicando por trimestre, con frecuencia de 4:

s.ts.q <- ts(s, start = 2001, frequency = 4)
plot(s.ts.q)

Podemos hacerle preguntas a una serie temporal:

start(s.ts.m) #♣ cuando empieza
## [1] 2001    1
end(s.ts.m) # cuando acaba
## [1] 2009    4
frequency(s.ts.m) # que frecuencia
## [1] 12

Trabajando con varias variables

Vamos a trabajar con varias variables usando precios:

prices <- read.csv("../DataSets/prices.csv")
head(prices)

Creamos el objeto temporal:

prices.ts <- ts(prices, start = c(1980,1), frequency = 12)
plot(prices.ts)

Vamos a juntar los dos gráficos en uno solo para ver más información:

plot(prices.ts, plot.type = "single", col = 1:2)
legend("topleft", colnames(prices.ts), col = 1:2, lty = 1)

Si se quieren añadir end para saber cuando acaba la serie temporal. La frecuencia es un elemento esencial, en el caso que se tomen medidas cada 10 minutos pondríamos frecuency = 6, con un start a la hora correspondiente.

La descomposición de una serie temporal

Una serie temporal se basa en encontrar una estacionalidado y una tendencia.

En este caso, parece que la amplitud de las fluctuaciones incrementa con el tiempo, también llamado seie temporal multiplicativa:

  • Para este caso dividiremos o aplicamos un logaritmo. Cuando tenemos la sospecha que tenemos una serie temporal multiplicativa, es decir el anterior valor multiplicado por un número, al aplicar el logaritmo esa propiedad se transforma en una suma:

    log (a*b) = log (a) + log (b)

Al hacer esta descomposición con log() podremos trabajar con dos factores para saber como afectan cada uno al precio.

Usaremos dos sistemas para hacer las series temporales:

  • stl() : Seasonal Decomposition of Time Series by Loess

  • decompose() : Classical Seasonal Decomposition by Moving Average. Va moviendo las medias y los promedios.

Primero pasaremos los valores a logaritmo para perder el factor multiplicativo y pasarlo a aditivo, descomponiendi así en suma el factor que dependa del tiempo y la tendencia, y el ruido:

#stl
#Seasonal Decomposition of Time Series by Loess
flour.l <- log(prices.ts[,1]) # seleccionamos la variable que  queremos estudiar
flour.stl<- stl(flour.l, s.window = "period") # usamos la función stl() para hacer la descomposición estacional, indicamos el s.window() para tomar el periodo original 
plot(flour.stl)

Podemos ver:

  • data : los valores en log(), negativos algunos, pero no se altera su forma

  • seasonal : factor estacional que se repite año a año

  • tendencia :

  • residuos: resto en el residuo, ruido, se puede hacer un histograma pero es ruta de la aleatoriedad.

Si sumasamos los residuos + tendencia + seasonal, sería igual a data.

Veamos cuáles son los resultados para el gas():

gas.l <- log(prices.ts[,2])
gas.stl <- stl(gas.l, s.window = "period")
plot(gas.stl)

Ahora trabajaremos con la descomposición clásica:

#decompose
#Classical Seasonal Decomposition by Moving Averages
flour.dec <- decompose(flour.l) #  se le aplicamos a los logaritmos de los precios
plot(flour.dec)

Vemos que hay ciertas diferencias respecto al anterior time series aplicado. Lo aplicamos sobre el dataset anterior.

gas.dec <- decompose(gas.l)
plot(gas.dec)

Podríamos analizar que ocurrió entre 2008 y 2009, tal y como nos indica el gráfico de “random”.

Ahora vamos a tomar los datos originales y ajustarlos, elliminando los datos estacionales para saber la tenencia de los datos. Tomamos la información de la gasolina:

gas.season.adjusted <- prices.ts[,2] - (gas.dec$seasonal)
plot(gas.season.adjusted)

El filtrado de series temporales para localizar tendencias

Usaremos la técnica de los promedios móviles, que toma un valor uallquiera de la serie temporal y se promedia con todos los anteriores.

s <- read.csv("../DataSets/ts-example.csv")

Definimos que el periodo de tiempo es 7, es decir que los datos son diarios, divididos en semanas:

n <- 7 # definimos un filtro por semanas, que es el que queremos analizar
weights <- rep(1/n, n) # crearemos un sistema de pesos de 7 en 7, que debem sumar entre todos 1

# Creamos el primer filtro 
s.fil.1 <- filter(s$sales, filter = weights, sides = 2) # iremos filtrando según pesos, y suvizaremos ya qu evamos haciendo el promedio de 3 anteriores, el actual, y los 3 posteriores. Indicamos esto con el parametro side = 2
plot(s$sales, type = "l")
lines(s.fil.1, col = "blue", lwd = 3)

Podemos hacer un filtrado unilateral, solo cogiendo los 6 anteriores, no 3 anteriores y 3 posteriores:

s.fil.2 <- filter(s$sales, filter = weights, sides = 1) # indicamos la dirección unidireccional de sides igual a 1

plot(s$sales, type = "l")
lines(s.fil.1, col = "blue", lwd = 3)
lines(s.fil.2, col = "red", lwd = 3)

Con estas dos téncincas podemos localizar las tendencias mucho mejor. Vamos aplicar estos casos sobre el dataset anteriore de los precios del gas:

n <- 12 # seleccionamos la ffrecuencia que es por meses
gas.f.1 <- filter(prices.ts[,2], filter = rep(1/n, n), sides = 2) # damos los pesos a los anteriores
gas.f.2 <- filter(prices.ts[,2], filter = rep(1/n,n), sides = 1) # damos distribución unilateral
plot(prices.ts[,2])
lines(gas.f.1, col = "blue", lwd = 3)
lines(gas.f.2, col = "red", lwd = 3)

Podemos ver la tendencia calculada de forma manual, sin usar decompose o stl.

Suavizado y predicción con el método de Holt-Winters

El paquete Holt-Winter lleva a cabo un suavizado exponencial en la presencia de tendencias y de la estacionalidad.El proceso a seguir para hacer una buena prediccións es:

  1. Eliminar el ruido de la muestra

  2. llevar a cavo un forecast, un apredicción de valores en el futuro

inf <- read.csv("../DataSets/infy-monthly.csv")
head(inf)

Creamos el objeto time series escogiendo el momento que empezo a recopilarse la información:

inf.ts <- ts(inf$Adj.Close, start=c(1999,3), frequency = 12)
inf.ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1999                      1.028697  0.985462  1.186729  1.383523  1.950755
## 2000  7.753204 12.899140  9.211965  9.689268  7.326615  8.463322  6.347494
## 2001  5.772664  4.393672  3.128103  3.510162  3.247497  3.107732  3.156978
## 2002  3.009222  2.633788  3.136603  2.971392  2.707535  2.456770  2.471164
## 2003  3.122863  2.970950  2.949317  1.975827  2.124855  2.588909  2.568642
## 2004  4.271798  4.062922  3.951715  3.903848  4.001034  4.577651  4.927434
## 2005  6.500268  7.505230  7.257188  5.827010  7.120370  7.638110  7.020751
## 2006  7.533632  6.997919  7.695732  7.773816  6.978148  7.645977  8.223350
## 2007 11.632120 10.882046 10.077827 10.498989  9.877274 10.103899  9.979263
## 2008  8.354819  7.854337  7.218644  8.816958  9.910752  8.906346  8.072274
## 2009  5.481445  4.994389  5.495892  6.358560  7.132483  7.590644  8.946654
## 2010 10.842394 11.884651 12.283591 12.519611 12.016239 12.590125 12.709910
## 2011 14.419378 14.204288 15.269077 13.880591 13.150146 13.990744 13.345149
## 2012 11.858818 12.438927 12.298751 10.200438  9.079036  9.859200  8.660170
## 2013 11.609458 11.875913 11.871508  9.191555  9.193759  9.178202 11.069996
## 2014 13.130185 13.822780 12.143964 12.038618 11.527576 12.052067 12.449975
## 2015 15.599781 16.803635 16.057522 14.180784 14.588176 14.510359 15.944217
## 2016 17.038885 16.001902 18.094896 17.885597 18.494471 16.981804 15.801635
##            Aug       Sep       Oct       Nov       Dec
## 1999  2.366856  3.412329  3.603229  4.915664  7.875515
## 2000  7.449423  6.302729  6.565342  5.723784  4.405611
## 2001  2.510091  1.616498  2.257170  2.734351  2.968998
## 2002  2.715882  2.600721  3.433240  4.002129  3.343520
## 2003  2.740914  3.284753  4.077592  4.015540  4.612671
## 2004  5.121977  5.561173  6.527009  7.086905  6.822131
## 2005  6.982283  7.326514  6.707096  7.122455  7.992251
## 2006  8.975842  9.552219 10.426787 10.735641 10.942215
## 2007  9.599007  9.735818 10.248864  8.506175  9.153975
## 2008  8.459593  6.826287  6.008607  5.190451  5.070749
## 2009  8.988236 10.081880  9.564166 10.646055 11.544195
## 2010 12.037428 14.145241 14.172558 14.087165 16.201830
## 2011 11.071626 10.953659 12.566576 11.132063 11.080306
## 2012  9.305633 10.620630  9.500367  9.788324  9.314873
## 2013 10.332439 10.720159 11.823147 12.108102 12.686385
## 2014 13.512834 13.737668 15.184336 15.984281 14.400500
## 2015 16.170509 16.425087 17.122826 15.868711 15.935307
## 2016 15.253432 15.176492 14.676380 14.072151 14.412293
plot(inf.ts)

Vamos a hacer el suavizado exponencial con un filtro diferente al que hemos usado anterioremente. El método que usa Holt-Winters es una convolución con un núcleo exponencial

inf.hw <- HoltWinters(inf.ts)
plot(inf.hw, col = "blue", col.predicted = "red")

Vemos en rojo como Holt-Winters ha suavizado los datos. Analizemos las variables que nos da Holt-Winters:

inf.hw$SSE # la suma de los cuadrados de los errores
## [1] 347.1345

Podemos ver los parámetros “alpha”,“beta” y “gamma” que se usan para llevar a cabo el filtro exponencial:

inf.hw$alpha
##     alpha 
## 0.5519187
inf.hw$beta
##       beta 
## 0.01042698
inf.hw$gamma # tiende a ir hacia 1
## gamma 
##     1

Consultando los valores ajustados(previsiones):

head(inf.hw$fitted)
##          xhat    level     trend      season
## [1,] 6.371099 5.038675 0.2702327  1.06219146
## [2,] 8.459056 6.876835 0.2865814  1.29563975
## [3,] 6.911948 7.842393 0.2936611 -1.22410629
## [4,] 8.684473 8.364917 0.2960474  0.02350837
## [5,] 6.968464 8.538907 0.2947748 -1.86521779
## [6,] 8.455786 8.490957 0.2912012 -0.32637179

Vamos a incorporar predicciones en base al modelo de Holt-Winters:

library(forecast)

infy.fore <- forecast(inf.hw, h=24) # podimos una predicción a dos años vista (24 meses)
plot(infy.fore)

Nos pinta una predicción donde se representa el 85% en azul oscuro y de 95% en azul claro.

Consultamos los valores del objeto creado para comprovar resultados:

head(infy.fore$lower) # valores más altos
##           80%      95%
## [1,] 13.82093 12.93681
## [2,] 13.44266 12.43035
## [3,] 13.87266 12.74443
## [4,] 12.02846 10.79309
## [5,] 11.46834 10.13248
## [6,] 10.25586  8.82474
head(infy.fore$upper) # valores más bajos
##           80%      95%
## [1,] 17.16121 18.04533
## [2,] 17.26724 18.27954
## [3,] 18.13522 19.26345
## [4,] 16.69579 17.93116
## [5,] 16.51534 17.85120
## [6,] 15.66275 17.09387

Creando un modelo autorregresivo integrado de media móvil

Los usaremos para series temporales univariantes. Usaremos el mismo dataset que antes

inf <- read.csv("../DataSets/infy-monthly.csv")
inf.ts <- ts(inf$Adj.Close, start = c(1999,3), frequency = 12)
inf.arima <- auto.arima(inf.ts)
summary(inf.arima)
## Series: inf.ts 
## ARIMA(2,1,1)(1,0,1)[12] 
## 
## Coefficients:
##           ar1      ar2     ma1     sar1    sma1
##       -0.8022  -0.0510  0.6271  -0.2878  0.3467
## s.e.   0.3108   0.1046  0.3023   0.9832  0.9679
## 
## sigma^2 estimated as 0.978:  log likelihood=-297.39
## AIC=606.78   AICc=607.19   BIC=626.95
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.06742221 0.9749815 0.7349418 0.4288637 9.862881 0.3035153
##                      ACF1
## Training set -0.004441285

El modelo arima se centra en encotrar 3 parámetros para encontrar la mejor predicción:

El módelo nos devuelve el criterio de acaice y bayesiano, que siempre deben ser lo más pequeño posible.

inf.fore <- forecast(inf.arima, h = 12)
plot(inf.fore, col = "red", # 
     fcol = "green")

Es un plot más útil a corto o medio plazo. En rojo está los datos reales y en verde la predicción.