Introduccion a Series de Tiempo
Eduardo Villarreal
8/5/2021
require(tidyverse)
require(gridExtra) #Libreria para hacer paneles de gráficas con ggplot2
require(forecast) #Libreria para análisis de series de tiempo
require(fpp2) #Datos de eries de tiempo1 Introducción
En esta sección estudiaremos los conceptos básicos de series de tiempo y utilizaremos la librería forecast de Rob Hyndman para modelar el comportamiento de las series de tiempos.
El forecasting es una tarea común en el ámbito de las ciencias, los negocios y la economía. Sin embargo, dependiendo de la aplicación, suele confudirse con procesos de negocios y en consecuencia es utilizado en la forma equivocada:
Predicción : se refiere a la tarea de predecir el futuro de la manera más exacata posible con la información disponible hasta el momento incluyendo información histórica cómo eventos futuros que puedan impactar la predicción.
Planeación: el objetivo es tomar las decisiones apropiadas para cumplir con una meta
Metas: es el resultado deseable como consecuencia de un proceso de planeación de recursos y la visión del futuro que se tiene de un negocio.
La tarea de hacer un pronóstico puede presentarse en 3 scenarios distitos:
Corto Plazo
Mediano Plazo
Largo Plazo
Y dependiendo del horizonte de planeación, la tarea de hacer un Forecast puede cambiar.
2 Forecast Package
El estándar para trabajar con series de tiempo dentro del forecast package es el objeto de clase ts el cual tiene propiedades distintas a los data frames. Para esta sección necesitamos instalar las librerias forecasty fpp2. Esta última contiene los datos necesarios.
2.1 Visualización de Series de tiempo
Una serie de tiempos es la combinación de 3 elementos:
Nivel
Tendencia
Estacionalidad
\(y_t = N_t + T_t + S_t +\epsilon_t\)
Usaremos los datos de venta de una medicina contra la diabetes para ejemplificar estos conceptos. El comando autoplot de la libreria forecast, utiliza como base ggplot para ayudar a visualizar la información. Los datos se encuentran en el data frame a10 de la libreria fpp2.
glimpse(a10)## Time-Series [1:204] from 1992 to 2008: 3.53 3.18 3.25 3.61 3.57 ...#Graficar la serie de tiempo
autoplot(a10) + theme_bw() + labs(title = "Venta de Medicamentos para Diabetes") + ylab("$Mio")
En este caso existen 2 patrones muy claros: existe una tendencia en las ventas y hay una estacionalidad.
Podemos visualizar el comportamiento estacional de las ventas con una gráfica de estacionalidad utilizando el comando ggseasonplot
ggseasonplot(ts_object)
ggseasonplot(a10, year.labels = TRUE, year.labels.left = TRUE) +
ylab("$Mio") +
ggtitle("Seasonal plot: Venta de medicinas para la Diabetes") +
theme_bw()
Una forma útil de esta gráfica es utilizar coordenadas polares para ver las variaciones año vs año (YoY):
ggseasonplot(a10, polar = TRUE) + ylab("$Mio") +
ggtitle("Seasonal plot: Venta de medicinas para la Diabetes")+
theme_gray()
2.2 Autocorrelación
Así como la correlación mide el grado de asociación entre 2 variables, la autocorrelación midel el grado de asociación entre un variable y sus resagos a lo largo del tiempo. Es decir, con la autocorrelación lo que nis interesa saber es el grado de asociación entre
\((y_t,y_{t−1}), (y_t,y_{t−2}), ...(y_t,y_{t−k})\)
El data frame beer dentro de la libreria fpp2 contiene la historia de venta de cerveza en Australia. Qué tan correlacionada se encuentra la serie a lo largo del tiempo?
# El siguiente codigo calcula el indice de correlacion de la serie para 12
# rezagos o Lags
# Paso 1 Convertir a data frame
temp = as.data.frame(beer)
temp$x = as.numeric(temp$x)
# Paso 2 Crear los 12 rezagos
require(quantmod)## Loading required package: quantmod## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last## Loading required package: TTRfor (i in 1:12) {
l = Lag(temp$x, i)
# Actualizar temp
temp[, (i + 1)] = l
}
# Paso 3 Eliminar NA´s
temp = na.omit(temp)
# Paso 4 Calcular atriz de Correlación
r = cor(temp)
r = data.frame(cor = r[, 1], name = row.names(r))
# Graficar Correlacion de y con y-k
ggplot(r[-1, ], aes(x = factor(name), y = cor, fill = cor)) + geom_bar(stat = "identity") +
scale_fill_gradient2(low = "blue4", mid = "yellow", high = "red") + theme_bw() +
labs(title = "Correlacion de y con y-k")
Tal cual se muestra la gráfica, podemos observar rezagos grandes, especialmente en el Lag(12) que representa una alta correlación con V13. Esto quiere decir que:
\(y_t=f(y_{t−1},y_{t−2},...,y_{t−12})\)
A esto se le conoce como Función de Autocorrelación y puede obtenerse mediante el siguiente código:
beer = ts(beer, frequency = 12)
ggAcf(beer, lag = 15)
2.3 Ruido Blanco (White Noise)
Una serie de tiempo cuya ACF es 0 se denomina White Noise. Este tipo de series tienen una propiedad importante conocida como Estacionariedad: Se dice que un proceso es Estacionario débil si y solo si la media y la varianza de dicho proceso se mantienen constantes a lo largo del tiempo; es decir, la distribución de probabilidad del proceso es la misma a lo largo del tiempo. Esto quiere decir que si un proceso es estacionario, el mejor estimado que podemos hacer es obtener \(E[y_t]=μE[y_t]=μ\)
# Simular un proceso de Ruido Blanco
set.seed(30)
y = ts(rnorm(100))
autoplot(y) + ggtitle("White noise")
# Calcular ACF
ggAcf(y) + labs(title = "Autocorrelacion de un proceso con Ruido Blanco")
Recuerdan las cartas de control de Shwart? Las cartas de contrl estan basadas en el supuesto de que la muestra de los datos sigue un proceso \(y_t\) que es estacionario y que además dicho proceso proviene de una distribución normal por lo que su mejor estimado es \(E[y_t]\) . Las cartas de control están diseñanadas para detectar movimientos en la media y varianza de un proceso y cuando esto sucede, se dice que el proceso no esta en control y requiere obaservación para identificar y eliminar causas de variación especial dentro del proceso.
Sin embargo, cuando observamos fenómenos como las ventas de cerveza, hablar de estacionariedad se vuelve complicado. Entoces la ACF ayuda a identificar patrones repetitivos en el tiempo para formular y proponer modelos de series de tiempos.
3 Operaciones con Series de tiempo
Muchas veces tenemos una serie de tiempos diaria y queremos hacer un summary por semana, por mes, quarto o año. Estos computos se pueden simplificar utilizando la libreria timetk cn la función `summarise_by_time`
Por ejemplo, vamos a obtener los precios de la acción de Apple, AAPL, y vamos a hacer un summary mensual de los precios de cierre:
require(timetk)## Loading required package: timetk## Registered S3 method overwritten by 'tune':
## method from
## required_pkgs.model_spec parsniprequire(tidyquant)## Loading required package: tidyquant## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union## Loading required package: PerformanceAnalytics##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legend## == Need to Learn tidyquant? ====================================================
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>#Descargar los datos de APPLE desde Yahoo Finance
getSymbols("AAPL", from = '2000-01-01',
to = "2021-05-01",warnings = FALSE,
auto.assign = TRUE)## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.## [1] "AAPL"#Graficar
AAPL = as.data.frame(AAPL)
AAPL$Date = row.names(AAPL)
AAPL$Date = as.Date(AAPL$Date, format = '%Y-%m-%d')
AAPL %>%
ggplot(aes(x = Date, y = AAPL.Close)) +
geom_line() +
labs(title = 'Precio de la Accion de Apple-AAPL') +
theme_bw()
#Agregar la serie de tiempos por MES
AAPL %>%
summarise_by_time(
Date, #Campo de la Fecha
.by = 'month', #Agregar por Mes
PriceClose = LAST(AAPL.Adjusted) #Obtener el ultimo precio del mes
) %>%
ggplot(aes(x = Date, y = PriceClose)) +
geom_line(color = 'blue', size = 1) +
theme_bw() +
labs(title = 'closing Price AAPL')
Podemos también obtener el promedio del precio por trimestre:
#Agregar la serie de tiempos por MES
AAPL %>%
summarise_by_time(
Date, #Campo de la Fecha
.by = 'quarter', #Agregar por Mes
PriceClose = AVERAGE(AAPL.Adjusted) #Obtener el ultimo precio del mes
) %>%
ggplot(aes(x = Date, y = PriceClose)) +
geom_line(color = 'blue', size = 1) +
theme_bw() +
labs(title = 'closing Price AAPL by Q')
4 Rolling Computations
Este tipo de computos se refiere a ventanas de tiempo bajo las cuales se calcula algún tipo de estadístico como la media la desviación estándar o la correlación entre dos series de tiempo. Un ejemplo es el cómputo de Media Móvil o MA(L):
\(MA(L) = \frac{1}{L}\sum_{t=1}^L X_t\)
Esto se puede generalizar para ventanas de tiempo con un ciclo for
for(i in L:T){
x_hat = mean(x[i:i+L, ])
}Para simplificar los computos vamos a utilizar la librería tidyquant que nos permitirá obtener datos de cualquier compañía listada en Yahoo Finance y hacer los cómputos para medias móviles, la volatilidad del precio y la correlación móvil entre el precio de cierre y el volúmen
4.1 Rolling Moving Average
Para este ejemplo vamos a utilizar datos del precio de la acción de Apple por medio de su ticker symbol AAPL. Para descargar los datos del precio de la acción, vamos a usar el comando getSymbols que toma como argumentos el ticker, la fecha de inicio y la fecha de fin.
require(tidyquant)
#Descargar los datos de APPLE desde Yahoo Finance
getSymbols("AAPL", from = '2000-01-01',
to = "2021-05-01",warnings = FALSE,
auto.assign = TRUE)## [1] "AAPL"#Graficar
AAPL = as.data.frame(AAPL)
AAPL$Date = row.names(AAPL)
AAPL$Date = as.Date(AAPL$Date, format = '%Y-%m-%d')
AAPL %>%
ggplot(aes(x = Date, y = AAPL.Close)) +
geom_line() +
labs(title = 'Precio de la Accion de Apple-AAPL') +
theme_bw()
Un procedimiento muy utilizado es usar dos medias móviles: una de un periodo corto (Fast MA) y otra de un periodo más largo (Slow MA) y observar el cruce entre ambas medias móviles para obtener señales de compra y/o venta.
En este ejemplo vamos a computar MA(L=28) y MA(L=80) :
#computar el promedio movil de 28 y 80 dias
AAPL = AAPL %>%tq_mutate(
# tq_mutate args
select = AAPL.Close,
mutate_fun = rollapply,
# rollapply args
width = 28,
align = "right",
FUN = mean,
# mean args
na.rm = TRUE,
# tq_mutate args
col_rename = "mean_28"
) %>%
tq_mutate(
# tq_mutate args
select = AAPL.Close,
mutate_fun = rollapply,
# rollapply args
width = 80,
align = "right",
FUN = mean,
# mean args
na.rm = TRUE,
# tq_mutate args
col_rename = "mean_80"
)
#Graficar los datos
AAPL %>%
filter(Date >= '2019-01-01') %>%
ggplot(aes(x = Date, y = AAPL.Close)) +
geom_line(color = 'blue') +
labs(title = 'Precio de la Accion de Apple-AAPL') +
theme_bw() +
geom_line(aes(y = mean_80), color = 'red') +
geom_line(aes(y = mean_28), color = 'orange')
4.2 Rolling Standard Deviation
Así como computamos la media móvil, podemos calcular la desviación estándar móvil:
\(Sd(L) = \frac{1}{L-1}\sqrt{\sum_{t=1}^L (X_t - \bar{X})^2}\)
AAPL = AAPL %>%tq_mutate(
# tq_mutate args
select = AAPL.Close,
mutate_fun = rollapply,
# rollapply args
width = 30,
align = "right",
FUN = sd,
# mean args
na.rm = TRUE,
# tq_mutate args
col_rename = "std_30"
)
#Graficar
AAPL %>%
filter(Date >= '2019-01-01') %>%
ggplot(aes(x = Date, y = std_30)) +
geom_line(color = 'blue') +
labs(title = 'STD de la Accion de Apple-AAPL') +
theme_bw()
Ahora ponemos todo en un solo panel con gridExtra
p1 = AAPL %>%
filter(Date >= '2019-01-01') %>%
ggplot(aes(x = Date, y = AAPL.Close)) +
geom_line(color = 'blue') +
labs(title = 'Precio de la Accion de Apple-AAPL') +
theme_bw() +
geom_line(aes(y = mean_80), color = 'red') +
geom_line(aes(y = mean_28), color = 'orange')
p2 = AAPL %>%
filter(Date >= '2019-01-01') %>%
ggplot(aes(x = Date, y = std_30)) +
geom_line(color = 'blue') +
labs(title = 'STD de la Accion de Apple-AAPL') +
theme_bw()
grid.arrange(p1, p2)
4.3 Rolling Correlation
Un supuesto importante de la correlación es que asumimos que en el tiempo se mantiene el tipo de relación aún cuando no es necesariamente cierto. Entonces, también podemos computar la correlación en ventanas de tiempo para 2 variables:
AAPL = AAPL %>%
tq_mutate_xy(
x = AAPL.Adjusted,
y = AAPL.Volume,
mutate_fun = runCor,
# runCor args
n = 252,
use = "pairwise.complete.obs",
# tq_mutate args
col_rename = "rolling_corr"
)
#Graficar
p3 = AAPL %>%
filter(Date >= '2019-01-01') %>%
ggplot(aes(x = Date, y = rolling_corr)) +
geom_line(color = 'blue') +
labs(title = 'Correlacion de la Accion de Apple-AAPL vs Volumen') +
theme_bw()
grid.arrange(p1, p2, p3)
En este caso, hemos computdo la correlación temporal para 1 año o 252 días. Es interesante ver cómo en el caso de la relación Precio de la Acción y volumen de la Acción hay mementos en el tiempo en donde la correlación es positiva!! Es decir, a mayor Precio mayor Volumen de transacciones y hay momentos en donde la correlación es negativa. Esto sugiere que la relación entre ambas variables no es estática como uno pensaría.
5 Descomposición de Series de Tiempo
Como vimos en la introducción, una serie de tiempos es la combinación de 3 elementos:
Nivel
Tendencia
Estacionalidad
\(y_t = N_t + T_t + S_t +\epsilon_t\)
5.1 Ajuste por Estacionalidad
Si removemos el componente estacional de los datos el resultado es una serie desestacionalida:
\(y_t - S_t = N_t + T_t + \epsilon_t\) si la series es aditiva y si la series es múltiplicativa, obtenemos:
\(\frac{y_t}{S_t} = N_t \times T_t \times \epsilon_t\)
Desestacionalizar una serie es útil cuando la variación estacional de la serie no es de interés. Por ejemplo, la tasa mensual de desempleo suele analizarse sin estacionalidad. Por ejemplo, una caída en la tasa de desempleo debido a los trabajos escolares.
Siguiendo con la serie de medicamentos para la diabetes a10, para ajustar por estacionalidad primero aplicamos la función decompose y posteriormente aplicamos la función seasadj:
#Descomponer la serie de tiempos
a10_decomp = decompose(a10)
#Ajustar por estacionalidad
a10_seasadj = seasadj(a10_decomp)
#Graficar
autoplot(a10_seasadj, color = 'red') +
geom_line(aes(y = a10), color = 'blue', lty = 2) +
theme_bw() +
labs(title = 'Serie desestacionalizada a10')
5.2 Ajuste por Tendencia
Si removemos el componente de tendencia de los datos el resultado es una serie sin tendencia pero con estacionalidad:
\(y_t - T_t = N_t + S_t + \epsilon_t\) si la series es aditiva y si la series es múltiplicativa, obtenemos:
\(\frac{y_t}{T_t} = N_t \times S_t \times \epsilon_t\)
Quitar la tendencia de la serie permite observar el nivel o valor promedio de la serie:
#Descomponer la serie de tiempos
a10_decomp = decompose(a10)
#Extraer la tendencia
trend = a10_decomp$trend
#Restar la tendencia a la serie original
a10_detrend = a10 - trend
#Graficar los resultados
autoplot(a10, color = 'red') +
geom_line(aes(y = a10_detrend), color = 'blue', lty = 2) +
theme_bw() +
labs(title = 'Serie sin tendencia a10')## Warning: Removed 12 row(s) containing missing values (geom_path).
5.3 Ajuste por Tendencia y Estacionalidad
Muchas veces queremos ver la serie eliminando la tendencia y la estacionalidad para poder observar el valor promedio de la serie:
#Descomponer la serie de tiempos
a10_decomp = decompose(a10)
#Extraer la tendencia y la estacionalidad
trend = a10_decomp$trend
seasonal = a10_decomp$seasonal
#Restar la tendencia a la serie original
a10_sesDetrend = a10 - trend - seasonal
#Graficar los resultados
autoplot(a10_sesDetrend) +
theme_bw() +
labs(title = 'Serie sin tendencia y estacionalidad a10')
La función decompose toma un objeto ts y arroja los componentes de la serie de forma individual:
a10_decomposed = decompose(a10)
autoplot(a10_decomposed)
6 Descomposición STL
La descomposición STL (Seasonal and Trend decomposition using Loess) tiene ventajas sobre otros métodos:
- Puede manejar diversos factores estacionales (mensual, trimestral, etc)
- El componente estacional puede cambiar en el tiempo
- El factor cíclico también puede cambiar en el tiempo
- Es un método robusto a los outliers
La descomposición STL puede hacerse con la función stl que toma como argumentos un vector de series de tiempo, t.window que se refiere al componente cíclico de la tendencia o el número de observaciones consecutivas para estimar el componente de tendencia y s.window que se refiere al factor de estacionalidad.
a10_stl = stl(a10, s.window="periodic", robust=TRUE, t.window = 13)
autoplot(a10_stl)