Clase 3

Series de tiempo en R: evaluando modelos basados en características

En la Clase 2 revisamos los aspectos teóricos básicos que nos permiten caracterizar a las series de tiempo.

En esta clase revisaremos las series de tiempo desde una perspectiva práctica, estableciendo los requerimientos necesarios para lograr manipularlas, agregarlas y pronosticarlas usando ecosistemas de R. Luego de esta clase podrás:

• Conocer los requerimientos técnicos (packages) para trabajar con series de tiempo en R
• Manipular y agregar series de tiempo
• Entender el proceso de evaluación de pronósticos a través del tiempo
• Definir indicadores de desempeño adecuados para evaluar modelos de pronóstico
• Entender la importancia de establecer una modelo ingenuo como pronóstico inicial
• Generar características para modelos de pronóstico a partir del índice temporal asociado a la serie
• Modelar series de tiempo utilizando Regresión Lineal

Objetos de R para Series de Tiempo

Desde la perspectiva programática, una serie de tiempo es una serie de observaciones cuyo principal atributo es poseer un índice temporal (timestamp) asociado a cada una de ellas. Dicho índice puede tomar la forma de un objeto de R tipo date, un objeto temporal u otro formato dependiendo de la frecuencia de la serie. La generación de dicho índice a partir de datos raw no es trivial y en general requerirá de un preprocesamiento o formateo de los datos a un formato de series de tiempo. Es por ello que conocer los ecosistemas que permitan al pracicante ser capaz de trabajar con períodos temporales y fechas se vuelve esencial.

base 🧓👴

Date

En ciertos lugares existen diferentes convenciones para trabajar con fechas (dates).

Show code

knitr::include_graphics("https://iso.mit.edu/wp-content/uploads/2020/01/am_dateformat.gif")

Figure 1: Fuente: Date Format in the United States

Sin embargo, existe un estándar global (ISO 8601 YYYY-MM-DD) que especifica la manera correcta para lidiar con fechas evitando toda confusión: ordenando los componentes de manera decreciente (años → meses → días). Cada valor posee una cantidad fija de digitos, por lo que es necesario rellenar con ceros algunos meses (i.e Septiembre: mes 9 → 09). La mejor manera de indicar a R que nos encontramos trabajando con fechas es declararlo de manera explícita utilizando el método as.Date sobre un string con una fecha en formato ISO.

Show code

2021-11-08 # No es una fecha

[1] 2002

Show code

str("2021-11-08") # String con una fecha

 chr "2021-11-08"

Show code

as.Date("2021-11-08") # Fecha

[1] "2021-11-08"

Show code

str(as.Date("2021-11-08")) # Estructura objeto Date

 Date[1:1], format: "2021-11-08"

Tras bambalinas, los objetos tipo Date son almacenados como los días transcurridos desde 1970-01-01, lo que significa que es posible llevar a cabo comparaciones matemáticas con ellas:

Show code

# ¿Es hoy mayor que ayer?
as.Date("2021-11-08") > as.Date("2021-11-07")

[1] TRUE

ademaś de restar unidades de tiempo,

Show code

# ¿Que fecha era hace un dia (una unidad de tiempo) atras?
as.Date("2021-11-08") - 1

[1] "2021-11-07"

y calcular diferencias entre fechas.

Show code

# ¿Cuantos dias han pasado desde la ultima clase?
as.Date("2021-11-08") - as.Date("2021-10-25")

Time difference of 14 days

Si disponemos de fechas, podemos hacer uso de la potencia de R para graficar con r-base:

Show code

x <- c(
  as.Date("2021-11-08"),
  as.Date("2021-10-08"),
  as.Date("2021-09-08"),
  as.Date("2021-08-08"),
  as.Date("2021-07-08")
)

set.seed(2021)
y <- rnorm(5)

plot(x, y, type = "b")

o bien con ggplot2

Show code

ggplot() +
  geom_line(aes(x = x, y = y)) +
  geom_point(aes(x = x, y = y))

Además, podemos obtener la fecha del sistema en el que nos encontremos trabajando con Sys.Date()

Show code

Sys.Date()

[1] "2022-01-19"

Time

Continuando la idea de ISO 8601, para incorporar un período en específico continuamos definiendo elementos de manera decreciente: HH:MM:SS, donde las horas poseen dos dígitos fijos (00 - 24) así como también los minutos (00 - 59) y segundos, pudiendo existir sin separador o con :.

Existen dos tipos de fechas en R:

• POSIXlt: lista con componentes nombrados
• POSIXct: segundos desde 1970-01-01 00:00:00, la coersión desde un string se lleva cabo mediante as.POSIXct()

Verifiquemos la estructura de un objeto POSIXct

Show code

lesson3_date <- as.POSIXct("2021-11-08 19:00:00")
str(lesson3_date)

 POSIXct[1:1], format: "2021-11-08 19:00:00"

Timezone

La ISO 8601 también permite especificar zonas, asumiendose zona local ante la ausencia de su definición:

Show code

# Coordinated Universal Time
as.POSIXct("2021-11-08 19:00:00", tz = "UTC")

[1] "2021-11-08 19:00:00 UTC"

La operatoria aritmética revisada con datos del tipo Date se extiende a los datos tipo POSIXct.

lubridate ⏰

A medida que los requerimientos de manipulación de datos han evolucionado a través del tiempo, también lo han hecho los paquetes con los cuales se abordan dichos desafíos. lubridate es un paquete de R que permite trabajar de manera fácil con tiempos y fechas¹.

Veamos como trabajar fechas a partir de distintos formatos:

Show code

lubridate::ymd("2021-11-08") # YYYY-MM-DD

[1] "2021-11-08"

Show code

lubridate::dmy("08/11/2021") # DD/MM/YYYY

[1] "2021-11-08"

Show code

lubridate::parse_date_time("Sep 11th, 2021", order = c("mdy"))

[1] "2021-09-11 UTC"

Show code

lubridate::parse_date_time("11th Sep 2021", order = c("dmy"))

[1] "2021-09-11 UTC"

Y obtener información de interés como:

Show code

# Anio
lubridate::year("2021-11-08")

[1] 2021

Show code

# Mes
lubridate::month("2021-11-08")

[1] 11

Show code

# Mes con etiqueta
lubridate::month("2021-11-08", label = TRUE) # factor!

[1] Nov
12 Levels: Jan < Feb < Mar < Apr < May < Jun < Jul < Aug < ... < Dec

Show code

# Dia
lubridate::day("2021-11-08")

[1] 8

Show code

# Dia de la semana
lubridate::wday("2021-11-08", label = TRUE) # factor!

[1] Mon
Levels: Sun < Mon < Tue < Wed < Thu < Fri < Sat

Show code

# Dia del anio
lubridate::yday("2021-11-08")

[1] 312

Show code

# Fecha-Hora-Minuto
lubridate::ymd_hm("2021-11-08 07:00pm")

[1] "2021-11-08 19:00:00 UTC"

Show code

# Hora
lubridate::hour("2021-11-08 07:00pm")

[1] 7

Show code

# Minuto
lubridate::minute("2021-11-08 07:00pm")

[1] 0

y la fecha en una zona en particular:

Show code

# Chile because of reasons
lubridate::with_tz("2021-11-08", "America/Santiago") 

[1] "2021-11-08 -03"

Además, podemos crear fechas con los métodos del tipo make_date*()

Show code

# Fecha
lubridate::make_date(year = 2021L, month = 11L, day = 8L)

[1] "2021-11-08"

Show code

# Tiempo
lubridate::make_datetime(year = 2021L, month = 11L, day = 8L, hour = 7L, min = 0L)

[1] "2021-11-08 07:00:00 UTC"

y obtener datos adicionales

Show code

# Trimestre
lubridate::quarter("2021-11-08 07:00pm")

[1] 4

Show code

# Semestre
lubridate::semester("2021-11-08 07:00pm")

[1] 2

Show code

# Anio bisiesto
lubridate::leap_year("2021-11-08 07:00pm")

[1] FALSE

Podemos restar fechas utilizando difftime()

Show code

difftime("2021-11-08", "2021-10-25", units = "weeks")

Time difference of 2 weeks

y obtener el ahora

Show code

lubridate::now()

[1] "2022-01-19 13:20:08 -03"

así como también la fecha de hoy.

Show code

lubridate::today()

[1] "2022-01-19"

En una semana más

Show code

lubridate::today() + lubridate::days(7)

[1] "2022-01-26"

En particular, cabe destacar dos definiciones en lubridate para intervalos temporales

• period: período o concepto humano de intervalo temporal. Una fecha (datetime) + un período de un día = mismo momento en la fecha siguiente.
• duration: duración o concepto cronometrado del tiempo. Una fecha + un período de un día = fecha + 86400 segundos.

Veamos la diferencia

Show code

lubridate::days(x = 7)

[1] "7d 0H 0M 0S"

Show code

lubridate::ddays(7)

[1] "604800s (~1 weeks)"

Lubridate CheatSheet

ts() 💊

Supongamos que tenemos los siguientes datos:

Año	Observación
2017	20
2018	30
2019	60
2020	40
2021	25

Podemos almacenar una serie de tiempo en un objeto tipo ts utilizando la función ts().

Show code

set.seed(2021)
y <- ts(rnorm(5), start=2017)

Es posible asociar observaciones que poseen una frecuencia mayor que la anual utilizando el argumento frequency

Show code

set.seed(2021)
y <- ts(rnorm(12*5), start=2017, frequency = 12)

En este caso, las frecuencias se definen como sigue

Dato	Frecuencia
Anual	1
Trimestral	4
Mensual	12
Semanal	522

¿Qué ocurre si la frecuencia de mis observaciones es mayor a la semanal?

En el caso de utilizar objetos ts, es necesario decidir que frecuencia es más representativa de la serie³ considerando efectos como la estacionalidad (en caso de que aplique).

xts() y zoo() 💊💊

El paquete xts (eXtensible time series) busca ofrecer un formato flexible y poderoso para trabajar con series de tiempo. Un objeto xts es un objeto tipo zoo (un índice + una matriz) extendido, donde cada fila corresponde a una observación en el tiempo.

Show code

x <- matrix(1:4, ncol = 2, nrow = 2)
idx <- as.Date(c("2021-10-01","2021-11-01"))
(x_xts <- xts::xts(x = x, order.by = idx))

           [,1] [,2]
2021-10-01    1    3
2021-11-01    2    4

Un aspecto relevante a notar es que el índice debe ser creciente en el tiempo, dejando las observaciones más recientes al final de la estructura de datos. Si se otorga un vector no ordenado, xts ordenará según el índice para asegurarse de que los datos estén apropiadamente ordenados. Si te fijas, un objeto xts parece poseer rownames, sin embargo, dichos valores corresponden al índice asociado a la serie (index). xts posee comportamientos especiales que incluyen: sus valores son matrices asociadas a fechas por cada observación, subconjuntos de los datos mantienen la forma matricial y además sus atributos son preservados.

En el caso de necesitar deconstruir un objeto xts, los métodos coredata() e index() nos permiten obtener sus componentes internos

Show code

zoo::coredata(x_xts, fmt = FALSE)

     [,1] [,2]
[1,]    1    3
[2,]    2    4

Show code

zoo::index(x_xts, fmt = FALSE)

[1] "2021-10-01" "2021-11-01"

tsibble() 💊💊💊

El paquete tsibble⁴ provee una infraestructura para datos temporales que adopta los principios de los datos tidy, siendo un objeto orientado a datos y modelos. Algunos aspectos importantes de los datos tipo tsibble son los siguientes:

1. Index es una variable con un orden inherente del pasado al presente
2. Key es un conjunto de variables que define unidades observacionales en el tiempo
3. Cada observación debe ser unicamente identificada por un index y su key
4. Cada unidad observacional debe ser medida bajo un intervalo común (si están regularmente espaciadas).

Show code

# Ejemplo de la documentacion de tsibble
weather <- nycflights13::weather %>% 
  select(origin, time_hour, temp, humid, precip)
weather_tsbl <- tsibble::as_tsibble(weather, key = origin, index = time_hour)
weather_tsbl %>% rmarkdown::paged_table()

Los objetos tsibble extienden a los data frames (tibble) introduciendo una estructura temporal que facilita el almacenamiento de múltiples series de tiempo en un sólo objeto, permitiendo aplicar funciones de dplyr como mutate(), filter(), select() y summarize() a los datos almacenados.

timetk() 💊💊💊💊🚀

Tal como hemos visto, existen muchos paquetes de R para trabajar con series de tiempo. Uno de los últimos ecosistemas desarrollados con gran éxito es modeltime⁵ que incluye paquete timetk. Podemos revisar una tabla comparativa de las virtudes de timetk en el siguiente enlace.

Revisemos algunas de las funcionalidades que otorga timetk para manipular series de tiempo⁶

Show code

tidyquant::FANG %>% rmarkdown::paged_table()

Podemos hacer facilmente gráficos de series de tiempo

Show code

tidyquant::FANG %>%
  group_by(symbol) %>%
  timetk::plot_time_series(date, adjusted, .facet_ncol = 2, .interactive = TRUE)

Resumir datos a través del tiempo

Show code

tidyquant::FANG %>%
  group_by(symbol) %>%
  timetk::summarise_by_time(
    date, .by = "quarter",
    volume = tidyquant::SUM(volume)
  ) %>%
  timetk::plot_time_series(date, volume, .facet_ncol = 2, .interactive = TRUE, .y_intercept = 0)

y suavizar (entre otras funcionalidades) series de tiempo:

Show code

tidyquant::FANG %>%
  group_by(symbol) %>%
  timetk::summarise_by_time(
    date, .by = "month",
    adjusted = tidyquant::FIRST(adjusted)
  ) %>%
  timetk::plot_time_series(date, adjusted, .facet_ncol = 2, .interactive = TRUE)

¿Entonces qué paquete debo utilizar? 🤔⁷.

Evaluación de pronósticos 📏

Set de entrenamiento y test 🔪

La precisión de un pronóstico sólo puede ser determinada considerando qué tan bien se desempeñan los pronósticos generados en datos nuevos que no han sido utilizados para ajustar el modelo (Hyndman and Athanasopoulos (2021)). Así, al igual que como se evalúan los modelos de Machine Learning, es una práctica común separar los datos disponibles en dos conjuntos, un conjunto de ajuste o entrenamiento (training set)⁸ y otro de evaluación (test set)⁹, donde los datos de entrenamiento son utilizados para estimar los parámetros de cualquier método de pronóstico y los datos de evaluación son utilizados para determinar su precisión. Dado que los datos de evaluación no son utilizados en el proceso de elaboración de los pronósticos, éstos nos permiten obtener un indicador confiable de qué tan bien puede pronosticar un modelo los datos nuevos.

Training Set

Test Set

Al momento de evaluar el desempeño de modelos de pronóstico es necesario distinguir entre un pronóstico (forecast) o valor predicho (predicted value) de \(y_t\), realizado en algún período previo, por ejemplo \(t-\tau\), y un valor ajustado (fitted value) de \(y_t\), que ha resultado de estimar los parámetros de un modelo de series de tiempo en datos históricos, generando residuales.

El método más común para evaluar el éxito de un pronóstico para predecir los valores reales es utilizar métricas de precisión o error. Aquí error no debe entenderse como una equivocación, sino como la parte impredecible de una observación. La selección de la métrica específica de error dependerá de los objetivos de pronóstico de que se tengan.

Error de pronóstico

El tamaño del set de evaluación depende de qué tanta información se posea y qué tan adelante se quiera pronosticar. Es ideal que el test set tenga al menos el largo del horizonte de pronóstico máximo requerido.

Algunas observaciones a considerar

• Un modelo que se ajusta bien a los datos de entrenamiento no necesariamente pronosticará bien.
• Siempre es posible obtener un ajuste perfecto usando un modelo con una cantidad de parámetros suficiente¹⁰
• Sobreajustarse (over-fit) un modelo a los datos es tan malo como fallar en identificar los patrones sistemáticos en ellos.

Validación Cruzada para Series de Tiempo (TSCV) 📝

Uno de los aspectos más relevantes del proceso de pronóstico planteado en la Clase 01 fue la validación de los modelos. En particular mencionamos la validación cruzada para series de tiempo¹¹.¹²

La validación cruzada de series de tiempo (a.k.a. backtesting) puede entenderse como una versión avanzada de la metodología de validación single out-of-sample utilizada para validar modelos de machine learning¹³ . Se basa en el uso de una ventana rodante (rolling window) para particionar la serie en múltiples pares de entrenamiento-evaluación. En general el proceso de validación cruzada involucra la creación de las particiones mencionadas, el entrenamiento o ajuste de un modelo con los datos de entrenamiento de cada partición, la evaluación de su desempeño por partición con los datos de evaluación de cada partición y la evaluación del modelo a nivel de precisión, escalabilidad y estabilidad basado en las métricas de desempeño obtenidas en los datos de evaluación por partición.

Como resultado del proceso anterior es posible generar un pronóstico final para revisar si se cumplen los criterios requeridos para su validez o bien aplicar ajustes adicionales y optimizar el modelo repitiendo el proceso de evaluación. Es usual considerar un modelo como estable si al examinar su desempeño en los datos de evaluación a través de las particiones, la distribución del error es estrecha, de este modo, los errores obtenidos en los verdaderos valores pronosticados deberían estar en el mismo rango de los errores de cada partición¹⁴.

Revisemos nuevamente el esquema de validación considerado en caret:

Show code

knitr::include_graphics("https://topepo.github.io/caret/splitting/Split_time-1.svg")

Figure 2: Fuente: Caret Package - Max Kuhn

De manera genérica¹⁵, los parámetros a considerar para llevar a cabo validación cruzada para series de tiempo son los siguientes:

• Largo de las particiones de entrenamiento o tamaño de ventana. Existen dos tipos
  • Ventana expansiva (expanding window): el tamaño de la partición de entrenamiento crece a través de las ventanas temporales¹⁶.
  • Ventana deslizante (sliding window): el tamaño de la partición de entrenamiento se mantiene fijo a través de las ventanas temporales¹⁷
• Largo de las particiones de evaluación u horizonte de pronóstico.
• Espacio¹⁸ entre particiones de entrenamiento y evaluación
• Cantidad de particiones

El método de la ventana expansiva es útil cuando la serie posee fuertes patrones estacionales y tendencia estable, pues las primeras observaciones podrían poseen información relevante que pueda ser utilizada por el modelo. La principal desventaja del método radica en que cada partición tendrá una cantidad diferente de datos de entrenamiento y los modelos aprenden mejor con mayor cantidad de datos¹⁹.

El sesgo anterior no está presente cuando se utiliza el método de la ventana deslizante, dado que todas las particiones poseen el mismo largo. Este método es recomendable cuando existen variaciones irregulares o cuando el poder predictivo está relacionado a la historia más reciente.

️🚨 Los beneficios de monitorear la estabilidad de un modelo conllevan un costo: la validación cruzada de series de tiempo²⁰ es un proceso computacionalmente costoso²¹.

Modelos benchmark (naive) 📐

En muchas organizaciones, la elaboración de pronósticos es llevada a cabo utilizando reglas y conocimiento de negocio. Este favorable escenario nos puede servir para establecer una situación base, benchmark o situación “sin proyecto.” Cualquier método de pronóstico que no provea de mejoras en rendimiento por sobre dichos escenarios base no será de valor para la organización y lo más lógico es que sea descartado. El caso contrario ocurre cuando no se dispone de un pronóstico preliminar o la tarea de pronóstico se encuentra siendo abordada por vez primera. En esta situación, antes de modelar y llevar a cabo esfuerzos de análisis y pronóstico de series de tiempo, es recomendable establecer una línea base sobre la cual construir soluciones más complejas²².

Algunos métodos de pronóstico pueden ser extremadamente simples y efectivos. Algunos métodos a considerar antes de realizar cualquier esfuerzo de pronóstico son los siguientes:

1. Pronosticar con el promedio histórico (mean method). Todos los valores futuros son iguales a la media de los datos históricos
2. Definir como pronóstico el valor asociado a la última observación. En este caso hablamos de un método ingenuo (naïve method)²³.
3. Definir como pronóstico el valor asociada a la última observación con estacionalidad similar. En este caso, hacemos uso del patrón estacional predominante en los datos, por ejemplo, asociando un pronóstico para el próximo mes equivalente al valor obtenido en la misma época del año pasado.

Veamos como podemos ajustar modelos ingenuos utilizando modeltime:

1. Cargamos librerías y datos

Show code

library(dplyr)
library(parsnip)
library(rsample)
library(timetk)
library(modeltime)
# Datos competencia M4
m750 <- m4_monthly %>% filter(id == "M750")
m750 %>% rmarkdown::paged_table()

2. Partición en conjuntos de entrenamiento y evaluación

Show code

# Particionamos los datos en proporcion 80/20
splits <- initial_time_split(m750, prop = 0.8)
splits

<Analysis/Assess/Total>
<244/62/306>

3. Definición de especificación de modelos

Show code

# Naive model ----
# Model Spec (parsnip)
model_spec_naive <- naive_reg(
  id = "id"
) %>%
  set_engine("naive")

model_spec_naive

Naive Forecast Model Specification (regression)

Main Arguments:
  id = id

Computational engine: naive 

Show code

# Seasonal Naive model ----
# Model Spec (parsnip)
model_spec_snaive <- naive_reg(
  id = "id",
  seasonal_period = 12
) %>%
  set_engine("snaive")

model_spec_snaive

Naive Forecast Model Specification (regression)

Main Arguments:
  id = id
  seasonal_period = 12

Computational engine: snaive 

4. Ajuste (entrenamiento) de modelos

Show code

# Ajuste ----
# Fit Spec
model_fit_snaive <- model_spec_snaive %>%
  fit(log(value) ~ date + id, data = training(splits))
model_fit_snaive

parsnip model object

Fit time:  40ms 
SNAIVE [12]
--------
Model: 
# A tibble: 12 × 3
   id    date       value
   <fct> <date>     <dbl>
 1 M750  2009-05-01  9.27
 2 M750  2009-06-01  9.27
 3 M750  2009-07-01  9.15
 4 M750  2009-08-01  9.19
 5 M750  2009-09-01  9.18
 6 M750  2009-10-01  9.25
 7 M750  2009-11-01  9.26
 8 M750  2009-12-01  9.27
 9 M750  2010-01-01  9.26
10 M750  2010-02-01  9.26
11 M750  2010-03-01  9.29
12 M750  2010-04-01  9.29

Show code

model_fit_naive <- model_spec_naive %>%
  fit(log(value) ~ date + id, data = training(splits))
model_fit_naive

parsnip model object

Fit time:  25ms 
NAIVE
--------
Model: 
# A tibble: 1 × 3
  id    date       value
  <fct> <date>     <dbl>
1 M750  2010-04-01  9.29

5. Establecimiento de tabla de modelos²⁴

Show code

# ---- MODELTIME TABLE ----
# Tibble con modelos ajustados
models_tbl <- modeltime_table(
  model_fit_naive,
  model_fit_snaive
)

6. Validación de modelo en datos de evaluación

Show code

# Validacion en test
calibration_tbl <- models_tbl %>%
  # Calibracion en datos de test para evaluacion de rendimiento
  modeltime_calibrate(new_data = testing(splits)) %>%
  # Pronostico en datos de test
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = m750
  )

7. Gráfico de pronóstico

Show code

calibration_tbl %>% 
# Grafico pronosticos
  plot_modeltime_forecast(.interactive = TRUE)

Métricas de desempeño 💯

MAE

En el error absoluto medio (Mean Absolute Error (MAE)), el error de pronóstico está en la misma escala que los datos originales. Debido a lo anterior, no puede ser utilizado para comparar desempeño entre series que involucren distintas unidades de medida. Un método de pronóstico que minimice el MAE conducirá a pronosticos de la mediana. Este método no es sensible a outliers.

\[MAE = \frac{1}{T}\sum_{t=1}^{T}\left|y_{t}-\hat{y}_{t}\right|\]

MSE

El error cuadrático medio (Mean Squared Error (MSE)) cuantifica la distancia al cuadrado promedio entre los valores reales y los pronosticados. El efecto de elevar al cuadrado previene que los valores negativos y positivos se cancelen entre sí, penalizando el error si es muy elevado. Un método de pronóstico que minimice el MSE conducirá a pronosticos de la media.

\[MSE = \frac{1}{T}\sum_{t=1}^{T}(y_t - \hat{y}_t)^2\]

RMSE

Similar a MSE, en las unidades de las observaciones.

\[RMSE = \sqrt{MSE} = \sqrt{\frac{1}{T}\sum_{t=1}^{T}(y_t - \hat{y}_t)^2}\]

MAPE

El error porcentual absoluto medio (Mean Absolute Percentage Error (MAPE)) es una de las métricas más fáciles de comparar y comunicar con interlocutores no-técnicos dado que representa un porcentaje. Los errores porcentuales poseen la ventaja de no tener unidades y son utilizados para comparar el desempeño de un modelo de pronóstico entre datasets. Sin embargo, poseen la desventaja de ser infintos o indefinidos si \(y_t = 0\) para cualquier periodo de interés \(t\), además de tener valores extremos cuando \(y_t\) es cercano a cero. Además, poseen la desventaja de penalizar errores negativos más que los positivos, lo que condujo a la elaboración de otra métrica porcentual llamada sMAPE (symmetric MAPE), usado en la competencia M3.

\[MAPE = \frac{1}{T}\sum _{t=1}^{T}\bigg(\frac{|y_{t}-\hat{y}_{t}|}{y_{t}}\bigg)\times100\]

SMAPE

\[SMAPE = {\frac{1}{T}\sum _{t=1}^{T}\bigg(\frac{2|y_{t}−\hat{y}_{t}|}{y_{t}+\hat{y}_{t}}\bigg)\times 100}\]

\(R^{2}\)

\[R^{2} = {1-{\frac {\sum _{t=1}^{T}\left(y_{t}-\hat{y}_{t}\right)}{\sum _{i=1}^{T}\left(y_{t}-\bar{y}_{t}\right)}}}\]

Obtengamos las métricas de rendimiento de los modelos ingenuos revisados anteriormente

Show code

# Metricas de rendimiento
models_tbl %>%
  # Calibracion en datos de test para evaluacion de rendimiento
  modeltime_calibrate(new_data = testing(splits)) %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(
    .interactive = TRUE
  )

Evaluación de desempeño en R 📚

Podemos encontrar las métricas revisadas anteriormente en los siguientes paquetes:

• Metrics
• yardstick

Regresión Lineal para pronóstico de series de tiempo 📈✨

La regresión lineal es uno de los métodos más utilizados para identificar y cuantificar la relación entre una variable dependiente²⁵ y una única²⁶ o múltiples²⁷ variables independientes²⁸.

Consideremos un modelo de regresión que permita una relación lineal entre el pronóstico \(y\) y un único predictor \(x\). El modelo de regresión lineal simple se puede escribir como

\[y_t = \beta_0 + \beta_{1}x_{t} + \epsilon_t\]

En este caso, serie de tiempo a predecir es \(y_t\), \(\beta_0\) es el intercepto²⁹, \(\beta_1\) es la pendiente que representa el cambio promedio de \(y\) ante el incremento unitario de \(x\), y \(\varepsilon_t\) es una serie de tiempo de residuos.

En el caso de la regresión lineal múltiple, se tiene más de una variable explicativa:

\[y_t = \beta_0 + \beta_{1}x_{1,t} + \beta_{2}x_{2,t} + \dots + \beta_{k}x_{k,t} + \epsilon_t\] Aquí, el modelo sólo incluye relaciones contemporáneas entre los regresores y la variable respuesta, por lo que hablamos de un modelo de regresión estático. Dichos modelos son apropiados cuando el valor esperado de la respuesta cambia inmediatamente cuando cambia el valor de una variable explicativa. Cada coeficiente de regresión \(\beta\) modela el cambio instantáneo en el valor esperado condicional de la variable \(y_t\) ante el cambio unitario de \(x_{k,t}\), dejando todas las demás preditores constantes.

Un modelo de regresión dinámico³⁰ incorpora relaciones entre valores actuales y rezagados de las variables independientes, lo que significa que la variable respuesta podría cambiar luego de un cambio en los valores de las variables explicativas.

\[\begin{aligned} y_t = \beta_0 & + \beta_{10}x_{1,t} + \beta_{11}x_{1,t-1} + ... + \beta_{1m}x_{1,t-m} \\ & + \beta_{20}x_{2,t} + \beta_{21}x_{2,t-1} + ... + \beta_{2m}x_{2,t-m} \\ & + \dots \\ & + \beta_{k0}x_{k,t} + \beta_{k1}x_{k,t-1} + ... + \beta_{km}x_{2,t-m} \\ & + \epsilon_t \\ \end{aligned}\]

Caso Aplicado

A continuación veremos cómo podemos analizar una serie de tiempo con R y ajustar un modelo de regresión lineal como aproximación para mejorar una aproximación ingenua. Para ello utilizaremos el dataset USgas que contiene la demanda de gas natural en US por estado y a nivel país.

Show code

library(USgas)
 
# The us_monthly dataset provides a monthly time series, 
# representing the demand for natural gas in the US between 2001 and 2020:
data("us_monthly")
head(us_monthly)

        date       y
1 2001-01-01 2676998
2 2001-02-01 2309464
3 2001-03-01 2246633
4 2001-04-01 1807170
5 2001-05-01 1522382
6 2001-06-01 1444378

• Particionamos los datos en proporcion 80/20

Show code

splits <- us_monthly %>% 
  timetk::tk_tbl() %>% 
  initial_time_split(prop = 0.8)

• Graficamos los datos de entrenamiento

Show code

training(splits) %>% 
  timetk::plot_time_series(
    .date_var = date,
    .value = y,
    .smooth = FALSE,
    .interactive = TRUE
  )

• Realizamos una descomposición en componentes de tendencia y estacionalidad

Show code

training(splits) %>% 
  plot_stl_diagnostics(
    .date_var = date, 
    .value = y,
    .feature_set = c("observed", "season", "trend", "remainder"),
    .interactive = TRUE)

• Ajustamos la tendencia mediante regresión lineal

Show code

training(splits) %>% 
  plot_time_series_regression(
    .date_var     = date,
    .formula      = y ~ as.numeric(date),
    .show_summary = TRUE,
    .facet_ncol   = 2,
    .interactive  = FALSE
  )

Call:
stats::lm(formula = .formula, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-536229 -301483 -160492  368938 1044090 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      816225.83  240717.41   3.391 0.000847 ***
as.numeric(date)     83.62      16.78   4.983  1.4e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 395300 on 191 degrees of freedom
Multiple R-squared:  0.115, Adjusted R-squared:  0.1104 
F-statistic: 24.83 on 1 and 191 DF,  p-value: 1.399e-06

• Ajustamos la estacionalidad mediante regresión lineal

Show code

training(splits) %>% 
  plot_time_series_regression(
    .date_var     = date,
    .formula      = y ~ lubridate::month(date, label = TRUE),
    .show_summary = TRUE,
    .facet_ncol   = 2,
    .interactive  = FALSE
  )

Call:
stats::lm(formula = .formula, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-591523 -135386  -34760  121129  447338 

Coefficients:
                                        Estimate Std. Error t value
(Intercept)                              2003361      13589 147.427
lubridate::month(date, label = TRUE).L   -502644      46894 -10.719
lubridate::month(date, label = TRUE).Q   1116715      46838  23.842
lubridate::month(date, label = TRUE).C     67621      46894   1.442
lubridate::month(date, label = TRUE)^4    165562      47000   3.523
lubridate::month(date, label = TRUE)^5    292944      47094   6.220
lubridate::month(date, label = TRUE)^6    -67106      47151  -1.423
lubridate::month(date, label = TRUE)^7   -186602      47178  -3.955
lubridate::month(date, label = TRUE)^8     71659      47186   1.519
lubridate::month(date, label = TRUE)^9     48622      47188   1.030
lubridate::month(date, label = TRUE)^10    74655      47189   1.582
lubridate::month(date, label = TRUE)^11    -7855      47189  -0.166
                                        Pr(>|t|)    
(Intercept)                              < 2e-16 ***
lubridate::month(date, label = TRUE).L   < 2e-16 ***
lubridate::month(date, label = TRUE).Q   < 2e-16 ***
lubridate::month(date, label = TRUE).C  0.151033    
lubridate::month(date, label = TRUE)^4  0.000541 ***
lubridate::month(date, label = TRUE)^5  3.35e-09 ***
lubridate::month(date, label = TRUE)^6  0.156399    
lubridate::month(date, label = TRUE)^7  0.000110 ***
lubridate::month(date, label = TRUE)^8  0.130597    
lubridate::month(date, label = TRUE)^9  0.304209    
lubridate::month(date, label = TRUE)^10 0.115385    
lubridate::month(date, label = TRUE)^11 0.867977    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 188800 on 181 degrees of freedom
Multiple R-squared:  0.8088,    Adjusted R-squared:  0.7972 
F-statistic: 69.62 on 11 and 181 DF,  p-value: < 2.2e-16

• Ajustamos la tendencia y estacionalidad mediante regresión lineal

Show code

training(splits) %>% 
  plot_time_series_regression(
    .date_var     = date,
    .formula      = y ~ as.numeric(date) + lubridate::month(date, label = TRUE),
    .show_summary = TRUE,
    .facet_ncol   = 2,
    .interactive  = FALSE
  )

Call:
stats::lm(formula = .formula, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-498244  -64561  -13618   74363  303095 

Coefficients:
                                          Estimate Std. Error t value
(Intercept)                              7.907e+05  7.096e+04  11.142
as.numeric(date)                         8.514e+01  4.947e+00  17.209
lubridate::month(date, label = TRUE).L  -5.265e+05  2.895e+04 -18.188
lubridate::month(date, label = TRUE).Q   1.109e+06  2.888e+04  38.391
lubridate::month(date, label = TRUE).C   7.485e+04  2.892e+04   2.589
lubridate::month(date, label = TRUE)^4   1.599e+05  2.898e+04   5.516
lubridate::month(date, label = TRUE)^5   2.970e+05  2.904e+04  10.228
lubridate::month(date, label = TRUE)^6  -6.961e+04  2.907e+04  -2.394
lubridate::month(date, label = TRUE)^7  -1.852e+05  2.909e+04  -6.369
lubridate::month(date, label = TRUE)^8   7.106e+04  2.909e+04   2.443
lubridate::month(date, label = TRUE)^9   4.889e+04  2.909e+04   1.680
lubridate::month(date, label = TRUE)^10  7.462e+04  2.909e+04   2.565
lubridate::month(date, label = TRUE)^11 -7.873e+03  2.909e+04  -0.271
                                        Pr(>|t|)    
(Intercept)                              < 2e-16 ***
as.numeric(date)                         < 2e-16 ***
lubridate::month(date, label = TRUE).L   < 2e-16 ***
lubridate::month(date, label = TRUE).Q   < 2e-16 ***
lubridate::month(date, label = TRUE).C    0.0104 *  
lubridate::month(date, label = TRUE)^4  1.19e-07 ***
lubridate::month(date, label = TRUE)^5   < 2e-16 ***
lubridate::month(date, label = TRUE)^6    0.0177 *  
lubridate::month(date, label = TRUE)^7  1.54e-09 ***
lubridate::month(date, label = TRUE)^8    0.0155 *  
lubridate::month(date, label = TRUE)^9    0.0946 .  
lubridate::month(date, label = TRUE)^10   0.0111 *  
lubridate::month(date, label = TRUE)^11   0.7870    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 116400 on 180 degrees of freedom
Multiple R-squared:  0.9277,    Adjusted R-squared:  0.9229 
F-statistic: 192.6 on 12 and 180 DF,  p-value: < 2.2e-16

• Ajustamos la tendencia de manera lineal y polinomial, además de la estacionalidad mediante regresión lineal

Show code

training(splits)  %>% 
  plot_time_series_regression(
    .date_var     = date,
    .formula      = y ~ as.numeric(date) + I(as.numeric(date)^2) + lubridate::month(date, label = TRUE),
    .show_summary = TRUE,
    .facet_ncol   = 2,
    .interactive  = FALSE
  )

Call:
stats::lm(formula = .formula, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-455876  -51547   -2126   56605  295319 

Coefficients:
                                          Estimate Std. Error t value
(Intercept)                              5.022e+06  5.756e+05   8.725
as.numeric(date)                        -5.175e+02  8.161e+01  -6.341
I(as.numeric(date)^2)                    2.116e-02  2.861e-03   7.395
lubridate::month(date, label = TRUE).L  -5.230e+05  2.541e+04 -20.584
lubridate::month(date, label = TRUE).Q   1.105e+06  2.535e+04  43.564
lubridate::month(date, label = TRUE).C   7.809e+04  2.538e+04   3.077
lubridate::month(date, label = TRUE)^4   1.573e+05  2.544e+04   6.182
lubridate::month(date, label = TRUE)^5   2.988e+05  2.549e+04  11.726
lubridate::month(date, label = TRUE)^6  -7.077e+04  2.552e+04  -2.774
lubridate::month(date, label = TRUE)^7  -1.846e+05  2.553e+04  -7.231
lubridate::month(date, label = TRUE)^8   7.075e+04  2.553e+04   2.771
lubridate::month(date, label = TRUE)^9   4.902e+04  2.553e+04   1.920
lubridate::month(date, label = TRUE)^10  7.458e+04  2.553e+04   2.921
lubridate::month(date, label = TRUE)^11 -7.863e+03  2.553e+04  -0.308
                                        Pr(>|t|)    
(Intercept)                             1.81e-15 ***
as.numeric(date)                        1.80e-09 ***
I(as.numeric(date)^2)                   5.25e-12 ***
lubridate::month(date, label = TRUE).L   < 2e-16 ***
lubridate::month(date, label = TRUE).Q   < 2e-16 ***
lubridate::month(date, label = TRUE).C   0.00242 ** 
lubridate::month(date, label = TRUE)^4  4.16e-09 ***
lubridate::month(date, label = TRUE)^5   < 2e-16 ***
lubridate::month(date, label = TRUE)^6   0.00613 ** 
lubridate::month(date, label = TRUE)^7  1.34e-11 ***
lubridate::month(date, label = TRUE)^8   0.00618 ** 
lubridate::month(date, label = TRUE)^9   0.05650 .  
lubridate::month(date, label = TRUE)^10  0.00394 ** 
lubridate::month(date, label = TRUE)^11  0.75849    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 102100 on 179 degrees of freedom
Multiple R-squared:  0.9446,    Adjusted R-squared:  0.9406 
F-statistic:   235 on 13 and 179 DF,  p-value: < 2.2e-16

• Ya tenemos una propuesta de modelo, podemos evaluar su capacidad predictiva out-sample

Show code

model_spec_lm <- linear_reg() %>%
  set_engine("lm")

model_fit_lm <- model_spec_lm %>%
  fit(
    y ~ as.numeric(date) + I(as.numeric(date) ^ 2) + lubridate::month(date, label = TRUE),
    data = training(splits)
  )

• Y, por supuesto, establecer un punto de comparación ingenuo

Show code

model_spec_snaive <- naive_reg(seasonal_period = 12) %>%
  set_engine("snaive")

model_fit_snaive <- model_spec_snaive %>%
  fit(y ~ date , data = training(splits))
model_fit_snaive

parsnip model object

Fit time:  19ms 
SNAIVE [12]
--------
Model: 
# A tibble: 12 × 2
   date         value
   <date>       <int>
 1 2016-02-01 2652260
 2 2016-03-01 2356298
 3 2016-04-01 2083848
 4 2016-05-01 1965799
 5 2016-06-01 2000656
 6 2016-07-01 2186616
 7 2016-08-01 2208375
 8 2016-09-01 1947752
 9 2016-10-01 1925203
10 2016-11-01 2159445
11 2016-12-01 2866273
12 2017-01-01 2913823

• Generamos la tabla de modelos³¹

Show code

models_tbl <- modeltime_table(
  model_fit_lm,
  model_fit_snaive
)
models_tbl

# Modeltime Table
# A tibble: 2 × 3
  .model_id .model   .model_desc
      <int> <list>   <chr>      
1         1 <fit[+]> LM         
2         2 <fit[+]> SNAIVE [12]

• Calibramos³² los modelos en el set de evaluación (test)

Show code

calibration_tbl <- models_tbl %>%
  # Calibracion en datos de test para evaluacion de rendimiento
  modeltime_calibrate(new_data = testing(splits)) %>%
  # Pronostico en datos de test
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = us_monthly
  )
calibration_tbl

# A tibble: 340 × 7
   .model_id .model_desc .key   .index      .value .conf_lo .conf_hi
       <int> <chr>       <fct>  <date>       <dbl>    <dbl>    <dbl>
 1        NA ACTUAL      actual 2001-01-01 2676998       NA       NA
 2        NA ACTUAL      actual 2001-02-01 2309464       NA       NA
 3        NA ACTUAL      actual 2001-03-01 2246633       NA       NA
 4        NA ACTUAL      actual 2001-04-01 1807170       NA       NA
 5        NA ACTUAL      actual 2001-05-01 1522382       NA       NA
 6        NA ACTUAL      actual 2001-06-01 1444378       NA       NA
 7        NA ACTUAL      actual 2001-07-01 1598071       NA       NA
 8        NA ACTUAL      actual 2001-08-01 1669178       NA       NA
 9        NA ACTUAL      actual 2001-09-01 1494128       NA       NA
10        NA ACTUAL      actual 2001-10-01 1649073       NA       NA
# … with 330 more rows

• Graficamos los pronósticos para comparar visualmente

Show code

calibration_tbl %>% 
  # Grafico pronosticos
  plot_modeltime_forecast(.interactive = TRUE)

• Revisamos el desempeño en datos de evaluación

Show code

models_tbl %>%
  # Calibracion en datos de test para evaluacion de rendimiento
  modeltime_calibrate(new_data = testing(splits)) %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(
    .interactive = TRUE
  )

Análisis de residuales 🚯

El análisis de residuales permite analizar qué tan bien el modelo captura e identifica los patrones de la serie. Además, podemos elaborar intervalos de confianza para los pronósticos generados a partir de la distribución de los residuales. Podemos definir a los residuales como la diferencia entre el valor real observado (\(y_i\)) y su correspondiente valor ajustado por el modelo (\(\hat{y_i}\)), para cada \(i = 1, \dots, T\).

\[\epsilon_i = y_i - \hat{y_i}\]

Revisemos los conceptos anteriores aplicados a los modelos ajustados para USgas:

• Generamos los residuales de entrenamiento

Show code

residuals_tbl <- models_tbl %>%
  modeltime_calibrate(new_data = training(splits)) %>%
  modeltime_residuals()

• Graficamos los residuales en el tiempo (entrenamiento)

Show code

residuals_tbl %>%
  plot_modeltime_residuals(
    .interactive = TRUE,
    .y_intercept = 0,
    .y_intercept_color = "blue"
  )

• Graficamos los función de autocorrelación de los residuales en el tiempo (entrenamiento)

Show code

residuals_tbl %>%
  plot_modeltime_residuals(
    .interactive = TRUE,
    .type = "acf"
  )

• Graficamos los estacionalidad (de existir) de los residuales en el tiempo (entrenamiento)

Show code

residuals_tbl %>%
  plot_modeltime_residuals(
    .interactive = TRUE,
    .type = "seasonality"
  )

• Obtenemos el error de pronóstico (test)

Show code

forecast_error_tbl <- models_tbl %>%
  modeltime_calibrate(new_data = testing(splits)) %>%
  modeltime_residuals()

• Graficamos el error de pronóstico en el tiempo (test)

Show code

forecast_error_tbl %>%
  plot_modeltime_residuals(
    .interactive = TRUE,
    .y_intercept = 0,
    .y_intercept_color = "blue"
  )

• Graficamos los función de autocorrelación del error en el tiempo (test)

Show code

forecast_error_tbl %>%
  plot_modeltime_residuals(
    .interactive = TRUE,
    .type = "acf"
  )

• Realizamos tests estadísticos a los residuales (entrenamiento)

Show code

models_tbl %>% 
  filter(.model_id == 1) %>% 
  modeltime_calibrate(new_data = training(splits)) %>%
  modeltime_residuals() %>%
  modeltime_residuals_test()

# A tibble: 1 × 6
  .model_id .model_desc shapiro_wilk box_pierce ljung_box
      <int> <chr>              <dbl>      <dbl>     <dbl>
1         1 LM              0.000931   5.27e-10  3.87e-10
# … with 1 more variable: durbin_watson <dbl>

• BONUS 🌀 : generación de características (variables) automática con timetk³³

Show code

training(splits) %>% 
  timetk::tk_augment_timeseries_signature(.date_var = date) %>% 
  timetk::tk_augment_fourier(date, .periods = 12, .K = 2) %>% 
  glimpse()

Rows: 193
Columns: 34
$ date          <date> 2001-01-01, 2001-02-01, 2001-03-01, 2001-04-0…
$ y             <int> 2676998, 2309464, 2246633, 1807170, 1522382, 1…
$ index.num     <dbl> 978307200, 980985600, 983404800, 986083200, 98…
$ diff          <dbl> NA, 2678400, 2419200, 2678400, 2592000, 267840…
$ year          <int> 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001…
$ year.iso      <int> 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001…
$ half          <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1…
$ quarter       <int> 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2…
$ month         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3…
$ month.xts     <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2,…
$ month.lbl     <ord> January, February, March, April, May, June, Ju…
$ day           <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ hour          <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ minute        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ second        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ hour12        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ am.pm         <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ wday          <int> 2, 5, 5, 1, 3, 6, 1, 4, 7, 2, 5, 7, 3, 6, 6, 2…
$ wday.xts      <int> 1, 4, 4, 0, 2, 5, 0, 3, 6, 1, 4, 6, 2, 5, 5, 1…
$ wday.lbl      <ord> Monday, Thursday, Thursday, Sunday, Tuesday, F…
$ mday          <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ qday          <int> 1, 32, 60, 1, 31, 62, 1, 32, 63, 1, 32, 62, 1,…
$ yday          <int> 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 3…
$ mweek         <int> 6, 5, 5, 5, 5, 5, 5, 5, 5, 6, 5, 5, 6, 5, 5, 6…
$ week          <int> 1, 5, 9, 13, 18, 22, 26, 31, 35, 40, 44, 48, 1…
$ week.iso      <int> 1, 5, 9, 13, 18, 22, 26, 31, 35, 40, 44, 48, 1…
$ week2         <int> 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1…
$ week3         <int> 1, 2, 0, 1, 0, 1, 2, 1, 2, 1, 2, 0, 1, 2, 0, 1…
$ week4         <int> 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 1, 1, 1, 1…
$ mday7         <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ date_sin12_K1 <dbl> 0.37877889, -0.13471146, -0.57126822, -0.90511…
$ date_cos12_K1 <dbl> -0.92548720, -0.99088487, -0.82076344, -0.4251…
$ date_sin12_K2 <dbl> -0.70111002, 0.26696710, 0.93775213, 0.7696512…
$ date_cos12_K2 <dbl> 0.7130531, 0.9637056, 0.3473053, -0.6384645, -…

Validacion Cruzada para series de tiempo

Show code

library(modeltime.resample) # TSCV
library(tidymodels)

# Especificacion Regresion Lineal-
model_spec_lm <- linear_reg() %>%
  set_engine("lm")

# Se extienden los datos originales con registros a ser pronosticados
us_monthly_full <- us_monthly %>%
  timetk::tk_tbl() %>%
  # Horizonte de pronostico h = 12
  timetk::future_frame(.date_var = date,
                       .length_out = "12 months",
                       .bind_data = TRUE)

# Revision de datos a pronosticar en el mismo objeto
us_monthly_full %>%
  tail(15)

# A tibble: 15 × 2
   date             y
   <date>       <int>
 1 2020-12-01 3158443
 2 2021-01-01 3286266
 3 2021-02-01 3036972
 4 2021-03-01      NA
 5 2021-04-01      NA
 6 2021-05-01      NA
 7 2021-06-01      NA
 8 2021-07-01      NA
 9 2021-08-01      NA
10 2021-09-01      NA
11 2021-10-01      NA
12 2021-11-01      NA
13 2021-12-01      NA
14 2022-01-01      NA
15 2022-02-01      NA

Show code

# Datos de entrenamiento
us_monthly_tbl <- us_monthly_full %>%
  filter(!is.na(y))

# Datos a pronosticar (fechas)
us_monthly_future_tbl <- us_monthly_full %>%
  filter(is.na(y))

# Cantidad de meses y anios para entrenar y validar
us_monthly_tbl %>%
  summarise(total_months = n_distinct(date)) %>%
  mutate(total_years = total_months / 12)

# A tibble: 1 × 2
  total_months total_years
         <int>       <dbl>
1          242        20.2

Show code

# Especificacion receta
# Forma funcional
recipe_spec <- recipes::recipe(y ~ date, data = us_monthly_full) %>%
  # Agrega Mes como variable categorica
  recipes::step_date(date, features = "month", ordinal = FALSE) %>%
  # Agrega tendencia lineal
  recipes::step_mutate(trend = as.numeric(date)) %>%
  # Agrega tendencia polinomial
  recipes::step_mutate(poly_trend = I(as.numeric(date)) ^ 2) %>%
  #recipes::step_normalize(date_num) %>%
  # Remueve la fecha para ajuste
  recipes::step_rm(date)

# Revisamos el eventual resultado de preprocesamiento
recipe_spec %>%
  prep() %>%
  juice() %>%
  tail(20)

# A tibble: 20 × 4
         y date_month trend poly_trend
     <int> <fct>      <dbl>   <I<dbl>>
 1 2492604 Jul        18444  340181136
 2 2404397 Aug        18475  341325625
 3 2174262 Sep        18506  342472036
 4 2323140 Oct        18536  343583296
 5 2440363 Nov        18567  344733489
 6 3158443 Dec        18597  345848409
 7 3286266 Jan        18628  347002384
 8 3036972 Feb        18659  348158281
 9      NA Mar        18687  349203969
10      NA Apr        18718  350363524
11      NA May        18748  351487504
12      NA Jun        18779  352650841
13      NA Jul        18809  353778481
14      NA Aug        18840  354945600
15      NA Sep        18871  356114641
16      NA Oct        18901  357247801
17      NA Nov        18932  358420624
18      NA Dec        18962  359557444
19      NA Jan        18993  360734049
20      NA Feb        19024  361912576

Show code

# Se establece un workflow que incorpora la receta de preprocesamiento y la
# la especificacion del modelo a ajustar. El workflow se ajusta para
# establecer la forma funcional
wflw_fit_lm <- workflow() %>%
  add_recipe(recipe_spec) %>%
  add_model(model_spec_lm) %>%
  fit(us_monthly_full)

# Esquema de validacion cruzada
us_monthly_tscv <- us_monthly_tbl %>%
  timetk::time_series_cv(
    date_var    = date,
    initial     = 12 * 15,
    # 15 anios
    assess      = "12 months",
    # Evaluacion
    skip        = "12 months",
    # Desplazamiento de ventana rodante
    cumulative  = TRUE # Acumulacion de datos historicos
  )

# Distancia temporal ascendente en particiones (slices)
# Grafico de esquema de validacion cruzada
us_monthly_tscv %>%
  timetk::tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(date, y,
                           .facet_ncol = 2, .interactive = TRUE)

Show code

# Tabla de modelos (Workflow modeltime)
(models_tbl <- modeltime_table(wflw_fit_lm))

# Modeltime Table
# A tibble: 1 × 3
  .model_id .model     .model_desc
      <int> <list>     <chr>      
1         1 <workflow> LM         

Show code

# Ajuste de modelo a distintas particiones
resamples_fitted <- models_tbl %>%
  modeltime.resample::modeltime_fit_resamples(resamples = us_monthly_tscv,
                                              control   = tune::control_resamples(verbose = TRUE))

── Fitting Resamples ────────────────────────────────────────────

5.569 sec elapsed

Show code

# Desempenio de modelos a traves de particiones (TSCV)
resamples_fitted %>%
  plot_modeltime_resamples(
    .point_size  = 3,
    .point_alpha = 0.8,
    .interactive = FALSE
  )

Show code

# Evaluacion de capacidad predictiva promedio
resamples_fitted %>%
  modeltime_resample_accuracy(summary_fns = mean) %>%
  table_modeltime_accuracy(.interactive = FALSE)

Accuracy Table
.model_id	.model_desc	.type	n	mae	mape	mase	smape	rmse	rsq
1	LM	Resamples	5	118230.4	4.9	0.47	4.79	148401.3	0.93

Show code

# Pronostico Real
resamples_fitted %>%
  modeltime_forecast(new_data = us_monthly_future_tbl,
                     actual_data = us_monthly_full) %>%
  plot_modeltime_forecast(.interactive = TRUE)

Show code

# Grafico comparativo particiones 
# Try-hard: podria ser un feature request para el paquete
slice_plot <- resamples_fitted$.resample_results %>%
  pluck(1) %>%
  mutate(
    training_tbl = map(splits, training),
    testing_tbl = map(splits, testing)
  ) %>%
  select(id, .predictions, training_tbl, testing_tbl) %>%
  mutate(assessment_tbl = pmap(list(training_tbl, testing_tbl, .predictions),
                               function(training_tbl, testing_tbl, .predictions) {
                                 testing_tbl %>%
                                   dplyr::bind_cols(.predictions %>% select(.pred)) %>%
                                   dplyr::bind_rows(training_tbl) %>%
                                   dplyr::mutate(.type = ifelse(is.na(.pred), "Training", "Testing")) %>%
                                   dplyr::arrange(date)
                                 
                               })) %>%
  select(id, assessment_tbl) %>%
  mutate(plot = map(assessment_tbl, function(x) {
    ggplot(data = x, aes(x = date, y = y, color = .type)) +
      geom_line() +
      geom_line(aes(x = date, y = .pred), color = "orange") +
      scale_color_manual(
        name = "Series",
        values = c(
          "Training" = "black",
          "Testing" = "lightblue",
          "Forecast" = "orange"
        )
      )
  }))

library(gridExtra)

grid.arrange(
  slice_plot$plot[[1]],
  slice_plot$plot[[2]],
  slice_plot$plot[[3]],
  slice_plot$plot[[4]],
  slice_plot$plot[[5]],
  ncol = 2
)

---

Show code

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Pop!_OS 21.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] gridExtra_2.3            vctrs_0.3.8             
 [3] rlang_0.4.12             yardstick_0.0.9         
 [5] workflowsets_0.1.0       workflows_0.2.4         
 [7] tune_0.1.6               recipes_0.1.17          
 [9] modeldata_0.1.1          infer_1.0.0             
[11] dials_0.0.10             scales_1.1.1            
[13] broom_0.7.10             tidymodels_0.1.4        
[15] modeltime.resample_0.2.0 USgas_0.1.1             
[17] modeltime_1.1.0          timetk_2.6.2            
[19] rsample_0.1.1            parsnip_0.1.7           
[21] xaringanExtra_0.5.5      forcats_0.5.1           
[23] stringr_1.4.0            dplyr_1.0.7             
[25] purrr_0.3.4              readr_2.1.1             
[27] tidyr_1.1.4              tibble_3.1.6            
[29] ggplot2_3.3.5            tidyverse_1.3.1         

loaded via a namespace (and not attached):
  [1] readxl_1.3.1               backports_1.4.0           
  [3] plyr_1.8.6                 lazyeval_0.2.2            
  [5] splines_4.1.2              crosstalk_1.2.0           
  [7] listenv_0.8.0              digest_0.6.29             
  [9] foreach_1.5.1              htmltools_0.5.2           
 [11] fansi_0.5.0                checkmate_2.0.0           
 [13] magrittr_2.0.1             memoise_2.0.1             
 [15] tzdb_0.2.0                 globals_0.14.0            
 [17] modelr_0.1.8               gower_0.2.2               
 [19] RcppParallel_5.1.4         tidyquant_1.0.3           
 [21] xts_0.12.1                 hardhat_0.1.6             
 [23] anytime_0.3.9              prettyunits_1.1.1         
 [25] colorspace_2.0-2           rvest_1.0.2               
 [27] haven_2.4.3                xfun_0.28                 
 [29] crayon_1.4.2               jsonlite_1.7.2            
 [31] progressr_0.9.0            survival_3.2-13           
 [33] zoo_1.8-9                  iterators_1.0.13          
 [35] glue_1.5.1                 gtable_0.3.0              
 [37] ipred_0.9-12               Quandl_2.11.0             
 [39] future.apply_1.8.1         quantmod_0.4.18           
 [41] DBI_1.1.1                  Rcpp_1.0.7                
 [43] viridisLite_0.4.0          GPfit_1.0-8               
 [45] lava_1.6.10                StanHeaders_2.21.0-7      
 [47] prodlim_2019.11.13         htmlwidgets_1.5.4         
 [49] httr_1.4.2                 ellipsis_0.3.2            
 [51] pkgconfig_2.0.3            farver_2.1.0              
 [53] nnet_7.3-17                sass_0.4.0.9000           
 [55] dbplyr_2.1.1               utf8_1.2.2                
 [57] tidyselect_1.1.1           labeling_0.4.2            
 [59] DiceDesign_1.9             munsell_0.5.0             
 [61] reactR_0.4.4               cellranger_1.1.0          
 [63] tools_4.1.2                cachem_1.0.6              
 [65] cli_3.1.0                  generics_0.1.1            
 [67] evaluate_0.14              fastmap_1.1.0             
 [69] yaml_2.2.1                 knitr_1.36                
 [71] fs_1.5.1                   future_1.23.0             
 [73] reactable_0.2.3            tictoc_1.0.1              
 [75] xml2_1.3.3                 compiler_4.1.2            
 [77] rstudioapi_0.13            plotly_4.10.0             
 [79] curl_4.3.2                 gt_0.3.1                  
 [81] reprex_2.0.1               lhs_1.1.3                 
 [83] bslib_0.3.1                stringi_1.7.6             
 [85] highr_0.9                  lattice_0.20-45           
 [87] Matrix_1.4-0               nycflights13_1.0.2        
 [89] pillar_1.6.4               lifecycle_1.0.1           
 [91] furrr_0.2.3                jquerylib_0.1.4           
 [93] data.table_1.14.2          R6_2.5.1                  
 [95] parallelly_1.29.0          codetools_0.2-18          
 [97] distill_1.3                MASS_7.3-55               
 [99] assertthat_0.2.1           withr_2.4.3               
[101] PerformanceAnalytics_2.0.4 parallel_4.1.2            
[103] hms_1.1.1                  tsibble_1.1.1             
[105] quadprog_1.5-8             grid_4.1.2                
[107] rpart_4.1-15               timeDate_3043.102         
[109] class_7.3-19               rmarkdown_2.11.3          
[111] downlit_0.4.0              TTR_0.24.2                
[113] pROC_1.18.0                lubridate_1.8.0           

---

1. Además de format parte del tidyverse↩︎

2. En un año hay \(365.25/7 = 52.18\) semanas en promedio, sin embargo, los objetos ts requieren de frecuencia dictadas por números enteros.↩︎

3. fpp2: ts objects↩︎

4. Hyndman et al: https://tsibble.tidyverts.org/↩︎

5. Matt Dancho’s modeltime↩︎

6. Puedes encontrar este ejemplo y otras viñetas en la documentación de timetk↩︎

7. A la fecha (2021-11-08), timetk se posiciona como uno de los mejores paquetes para manipulación y modelamiento de series de tiempo. Presenta la mayor cantidad de funcionalidades y wrappers en comparación a otros paquetes. El ecosistema de Matt Dancho (modeltime) está siendo extendido por otros colaboradores que han agregado más modelos como por ejemplo: garchmodels, bayesmodels y boostime. Sin embargo, esta imagen plantea bastantes desafíos en relación a la mantención de soluciones de este tipo. ¿Recomendación? ¡Aprender todos los formatos y saber moverse de uno a otro!📖↩︎

8. a.k.a. in-sample data↩︎

9. a.k.a. out-sample data, hold-out set↩︎

10. Por ejemplo, redes neuronales↩︎

11. fpp3: Time series cross-validation↩︎

12. caret: Data Splitting for Time Series↩︎

13. Un muestreo aleatorio de una serie de tiempo probablemente no sea la mejor idea para muestrear datos que presentan dependencias temporales.¡Ya sabemos que el orden temporal sí importa!↩︎

14. Asumiendo que no hay comportamientos anormales que impacten las métricas de desempeño consideradas.↩︎

15. Independiente del paquete utilizado↩︎

16. Véase fila inferior de la Figura 2↩︎

17. Véase fila superior de la Figura 2↩︎

18. de ser requerido↩︎

19. No es sorprendente notar que el rendimiento del modelo en las últimas particiones sea mejor que el de las primeras.↩︎

20. Y como veremos más adelante, la optimización de hiper-parámetros de modelo…↩︎

21. Pero por otro lado, sabemos que la validación cruzada es un proceso vergonzosamente paralelizable y actualmente tenemos a disposición alternativas de cómputo en la nube💅 ….↩︎

22. Lamentablemente, este no es un curso de problem solving 🙁. ¡La práctica le enseñará a no matar moscas con tanques! 🐷↩︎

23. Funciona sorprendentemente bien en series de tiempo financieras. También son llamadas caminatas aleatorias (random walk), como veremos más adelante.↩︎

24. Requerida por ecosistema modeltime↩︎

25. output, endógena,variable explicada,pronóstico, “\(y\)”↩︎

26. Regresión Lineal Univariada↩︎

27. Regresión Lineal Multivariada↩︎

28. input, exógena,variable explicativa, predictora, “\(x\)”↩︎

29. Valor constante que representa la media del valor esperado de \(y_t\) cuando \(x_t = 0\)↩︎

30. El término tambien es utilizado en la literatura para describir de manera general a los modelos de regresión con errores autocorrelacionados en el tiempo↩︎

31. Requerimiento de ecosistema modeltime↩︎

32. Terminología modeltime para ajuste y pronóstico de valores↩︎

33. ¡Próximas clases!↩︎