Series de Tiempo

Universidad Central

Dirección de Educación Continua.

Series de tiempo

Universidad Central

¿Qué se puede pronosticar?

Hay muchos fenómenos o variables que se pueden pronosticar; p. ej. los insumos que va a requerir una planta el siguiente mes o semana, las ventas de cierta sucursal en el mediano o largo plazo, la demanda de energía para una región del país en las siguientes 24 horas, etc.

Así, podemos necesitar pronósticos para muy corto plazo (minutos, horas,…), diarios, semanales, mensuales, hasta de largo plazo (para los siguientes 5, 10 años).

Es importante saber que existen variables que se pueden pronosticar con gran exactitud (p. ej., la puesta del sol el día de mañana, o el siguiente eclipse lunar), y otras que, por mas que intentemos, nuestro pronostico sera un mero experimento aleatorio (pronosticar los números ganadores de la lotería).

Entonces, que tan predecible es una variable o fenómeno depende de varias cosas:

1. Que tan bien conocemos los factores o elementos que influyen en él.
2. Que tantos datos tenemos disponibles.
3. Si el pronostico que vamos a realizar puede influir en el resultado de la variable a predecir.

Un caso de un tipo de series que son difíciles de pronosticar, son los tipos de cambio. Estos solo cumplen una de las tres características: tenemos muchos datos disponibles de ellos. Pero, desafortunadamente, hasta el momento no conocemos bien todo lo que afecta al tipo de cambio. Aparte, si se realiza algún pronostico de que, por ejemplo, el tipo de cambio va a aumentar, este puede influir en el valor futuro del mismo, logrando una “profecía autocumplida”. Esto es:

• Debido a las leyes de la oferta y la demanda, y a las nuevas expectativas de la gente a que va a subir el precio del dólar, la gente precavida comenzara a comprar dólares mientras están “baratos”.
• Esta nueva demanda ejercerá presión en el tipo de cambio y comenzará a aumentar.
• Al ver esto las demás personas, verán que se esta cumpliendo el pronostico y querrán adquirir mas dólares, antes de que suba tanto el precio.
• Esto podría provocar un efecto de bola de nieve y presionar mas el tipo de cambio.

Un buen pronóstico captura los patrones genuinos de la variable y las relaciones existentes en los datos históricos, pero no replica eventos pasados que no volverán a suceder. Otra cualidad de un buen modelo de pronóstico es aquel que captura la forma en que las cosas están cambiando.

Pronósticos, planeación y metas

En los negocios, realizar pronósticos es una tarea estadística muy común; p. ej.

• Las plantas ensambladoras requieren realizar pronósticos de la demanda, para agendar la cantidad de obreros a requerir por turno.
• Una sucursal de carnicería debe prever cuantos kilogramos de carne se van a demandar al día siguiente, para no quedarse sin existencias, ni con demasiados sobrantes.
• Una empresa de outsourcing de call center debe pronosticar la cantidad de llamadas que recibirá, para habilitar ejecutivos de atención a clientes suficientes.
• etc.

Sin embargo, los pronósticos de negocios a menudo se realizan muy deficientemente y son confundidos con la planeación y las metas de la empresa.

• Los pronósticos se tratan de predecir el futuro de la manera mas certera posible, dada la información disponible, incluyendo datos históricos y conocimiento de eventos futuros que puedan impactar el pronóstico.
• Los objetivos es lo que la empresa quiere que suceda (p. ej. crecer las ventas el 80%, YOY). Los objetivos o metas deberían estar ligados al pronóstico, pero no siempre suceden. Seguido las metas se establecen sin ningún plan para como lograrlo o sin un pronóstico de si son asequibles o no.
• La planeación es la respuesta a los pronósticos y metas. Involucra llevar a cabo las acciones necesarias para que los pronósticos igualen a las metas.

Depende de la aplicación específica, pero las empresas necesitan pronósticos de corto, mediano y largo plazo.

• Los pronósticos de corto plazo se requieren para la programación de personal, producción y transporte, principalmente. Para ello, pronósticos de la demanda son requeridos de igual manera.
• Los pronósticos de mediano plazo se utilizan para determinar los requerimientos futuros de recursos, para poder adquirir materia prima, contratación de personal, compra de maquinaria y equipo, etc.
• Los pronósticos de largo plazo son utilizados en la planeación estratégica. Tales decisiones deben tomar en cuenta oportunidades de mercado y factores internos y externos de la empresa.

Los pasos para llevar a cabo un pronóstico

1. La definición del problema. Se necesita tener un buen entendimiento del modo en que se usará el pronóstico, quién lo va a consumir y cómo encaja dentro de la organización.
2. Recopilación de información. Por lo general se requieren dos tipos de información: datos estadísticos y el expertise acumulado de la gente que recopila la información y usa los pronósticos.
3. Análisis exploratorio de los datos. Siempre se debe comenzar graficando los datos. ¿Se muestran patrones relevantes?, ¿existe una tendencia significativa?, ¿la estacionalidad es importante?, ¿existe evidencia de ciclos de negocio?, ¿hay datos outliers que se requieran explicar?
4. Selección y estimación de modelos. La selección del modelo dependerá de la disponibilidad de la información y del uso que se le dará al pronóstico. Es común comparar entre dos o tres modelos potenciales.
5. Uso y evaluación del modelo de pronóstico. Una vez seleccionado el modelo y estimado sus parámetros, se realiza el pronóstico. La evaluación del desempeño del modelo solo se puede evaluar hasta que se obtiene la información real del periodo pronosticado.

El análisis de series de tiempo

Para el análisis de series de tiempo, lo primero que se debe tener es la serie limpia. Una vez teniendo los datos listos para trabajar, la primera parte del análisis es visual: debemos graficar la o las serie(s).

Los gráficos nos permiten observar patrones, relaciones entre variables, cambios en el tiempo, valores inusuales, etc.

En R, podemos trabajar con objetos de series de tiempo, especialmente con la clase tsibble, que es un data frame consciente del tiempo.

Preparando los datos

Para la primera parte del entrenamiento, estaremos utilizando los datos que contienen la producción semanal de gasolina en los Estados Unidos desde el 3 de febrero de 1991 hasta el 15 de enero de 2017, en millones de barriles.

Cargando los datos

# Cargar los datos desde un archivo CSV
data_url <- "https://raw.githubusercontent.com/datacamp/time-series-analysis-in-python-live-training/master/data/US_gasoline_production.csv"
gasoline_prod_raw <- read_csv(data_url)

## Rows: 1355 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (1): value
## date (1): date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Imprimir información del DataFrame
glimpse(gasoline_prod_raw)

## Rows: 1,355
## Columns: 2
## $ date  <date> 1991-02-03, 1991-02-10, 1991-02-17, 1991-02-24, 1991-03-03, 199…
## $ value <dbl> 6.621, 6.433, 6.582, 7.224, 6.875, 6.947, 7.328, 6.777, 7.503, 6…

# Imprimir el DataFrame
head(gasoline_prod_raw)

## # A tibble: 6 × 2
##   date       value
##   <date>     <dbl>
## 1 1991-02-03  6.62
## 2 1991-02-10  6.43
## 3 1991-02-17  6.58
## 4 1991-02-24  7.22
## 5 1991-03-03  6.88
## 6 1991-03-10  6.95

Creación de un objeto de serie temporal (tsibble)

Los datos contienen fechas en la columna date, pero la columna es de tipo character. Para trabajar con series de tiempo en R de una manera moderna y ordenada, convertiremos este data frame en un objeto tsibble, especificando cuál es la columna de índice (la fecha) y la frecuencia de los datos.

# Convertir la columna de fecha a formato Date y luego a tsibble
gasoline_prod_weekly <- gasoline_prod_raw %>%
  mutate(date = as_date(date)) %>%
  as_tsibble(index = date)

# Imprimir un resumen del tsibble
glimpse(gasoline_prod_weekly)

## Rows: 1,355
## Columns: 2
## $ date  <date> 1991-02-03, 1991-02-10, 1991-02-17, 1991-02-24, 1991-03-03, 199…
## $ value <dbl> 6.621, 6.433, 6.582, 7.224, 6.875, 6.947, 7.328, 6.777, 7.503, 6…

# Imprimir el tsibble
head(gasoline_prod_weekly)

## # A tsibble: 6 x 2 [7D]
##   date       value
##   <date>     <dbl>
## 1 1991-02-03  6.62
## 2 1991-02-10  6.43
## 3 1991-02-17  6.58
## 4 1991-02-24  7.22
## 5 1991-03-03  6.88
## 6 1991-03-10  6.95

Análisis exploratorio

Graficando la serie de tiempo

# Graficar la serie usando autoplot para tsibbles
gasoline_prod_weekly %>%
  autoplot(value) +
  labs(
    title = "Producción Semanal de Gasolina en EE. UU.",
    x = "Fecha",
    y = "Millones de Barriles"
  )

El gráfico parece muy ruidoso. ¿Realmente necesitamos datos semanales?

Re-muestreo

El re-muestreo implica cambiar la frecuencia de los datos:

Upsampling: de una frecuencia más baja a una más alta (diaria → por hora). Downsampling: de una frecuencia más alta a una más baja (diaria → mensual).

Los datos de frecuencia más baja suelen agregarse: el valor para el mes X es el promedio/suma/máximo/… de todos los valores semanales en el mes X.

En R con tsibble, usamos index_by() para cambiar la frecuencia y summarise() para agregar.

# Re-muestrear los datos a valores promedio mensuales
gasoline_prod_monthly <- gasoline_prod_weekly %>%
  index_by(year_month = ~ yearmonth(.)) %>%
  summarise(value = mean(value))

# Imprimir el encabezado de los datos re-muestreados
head(gasoline_prod_monthly)

## # A tsibble: 6 x 2 [1M]
##   year_month value
##        <mth> <dbl>
## 1  1991 feb.  6.72
## 2  1991 mar.  7.09
## 3  1991 abr.  6.97
## 4  1991 may.  7.19
## 5  1991 jun.  7.62
## 6  1991 jul.  7.45

El tsibble re-muestreado ahora contiene el mes como índice, y el valor correspondiente es el promedio de los valores de ese mes.

# Graficar los datos re-muestreados
gasoline_prod_monthly %>%
  autoplot(value) +
  labs(
    title = "Producción Mensual Promedio de Gasolina en EE. UU.",
    x = "Fecha",
    y = "Millones de Barriles"
  )

Parece que hay algún tipo de patrón y también una tendencia, pero ambos son difíciles de ver.

Restringe los datos a un rango de tiempo para “acercar” la trama:

# Graficar los datos de 1999 a 2002
gasoline_prod_monthly %>%
  filter_index("1999" ~ "2002") %>%
  autoplot(value) +
  labs(
    title = "Producción Mensual Promedio (1999-2002)",
    x = "Fecha",
    y = "Millones de Barriles"
  )

## Warning: There were 2 warnings in `filter()`.
## The first warning was:
## ℹ In argument: `time_in(year_month, ...)`.
## Caused by warning:
## ! `yearmonth()` may yield unexpected results.
## ℹ Please use arg `format` to supply formats.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

Suavizado: ventanas deslizantes

Incluso después de re-muestrear, las tendencias y patrones podrían ser difíciles de ver.

La agregación deslizante (promedio móvil) crea una nueva serie temporal: el valor para cada marca de tiempo es una agregación de un número fijo de valores anteriores.

Ventana más grande = menos detalle = más suavidad

En R, podemos usar la función rollmean() del paquete zoo junto con mutate() de dplyr.

# Calcular el promedio móvil con una ventana de 2 meses
gasoline_prod_monthly %>%
  mutate(rolling_2m = rollmean(value, k = 2, fill = NA, align = "right")) %>%
  head()

## # A tsibble: 6 x 3 [1M]
##   year_month value rolling_2m
##        <mth> <dbl>      <dbl>
## 1  1991 feb.  6.72      NA   
## 2  1991 mar.  7.09       6.90
## 3  1991 abr.  6.97       7.03
## 4  1991 may.  7.19       7.08
## 5  1991 jun.  7.62       7.41
## 6  1991 jul.  7.45       7.54

Ejercicio

1. Graficar los datos remuestreados con un promedio móvil de 6 meses.

gasoline_prod_monthly %>%
  mutate(rolling_6m = rollmean(value, k = 6, fill = NA, align = "right")) %>%
  autoplot(value) +
  geom_line(aes(y = rolling_6m), color = "blue", size = 1) +
  labs(
    title = "Producción de Gasolina con Promedio Móvil de 6 Meses",
    x = "Fecha", y = "Millones de Barriles"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).

2. Graficar los datos remuestreados con un promedio móvil de 12 meses.

gasoline_prod_monthly %>%
  mutate(rolling_12m = rollmean(value, k = 12, fill = NA, align = "right")) %>%
  autoplot(value) +
  geom_line(aes(y = rolling_12m), color = "red", size = 1) +
  labs(
    title = "Producción de Gasolina con Promedio Móvil de 12 Meses",
    x = "Fecha", y = "Millones de Barriles"
  )

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_line()`).

Serie temporal = Tendencia + Estacionalidad + Residuos

• Tendencia: dirección general de la serie (por ejemplo, aumento)
• Estacionalidad: patrones repetitivos en los datos (por ejemplo, ventas durante las vacaciones…)
• Residuos: el resto

La función STL() del paquete feasts toma un tsibble y devuelve una descomposición de la serie en sus componentes.

# Descomponer la serie temporal usando STL
gasoline_prod_monthly %>%
  model(STL(value)) %>%
  components() %>%
  autoplot()

Nota: La función autoplot de feasts muestra la serie original y sus componentes de tendencia, estacionalidad y residuos.

No confíes solo en tus ojos

Una serie que muestra una tendencia o estacionalidad claramente visible probablemente realmente tiene esas características. ¡Pero si no puedes detectar la tendencia o la estacionalidad a simple vista, no significa que no estén presentes!

Estos datos se encuentran en una curva senoidal perturbada. ¿Puedes decirlo?..

La autocorrelación nos ayudará a identificar patrones ocultos en la serie.

Autocorrelación

Para hacer predicciones, necesitamos saber cómo los valores futuros de una serie temporal dependen de los valores pasados.

La correlación de dos variables aleatorias es una medida de la dependencia de una variable sobre la otra. Para las series temporales, la autocorrelación, como el nombre podría sugerir, medirá la autocorrelación, es decir, la dependencia de una serie temporal sobre sí misma.

Retardos (Lags)

Una serie temporal retardada es simplemente la misma serie temporal desplazada por un cierto número de días/semanas/meses, llamado el retardo (lag), en el pasado:

¡Observa los NA! Cuanto mayor es el retardo, más corta será la serie retrasada.

Función de autocorrelación

Función del retardo \(\mathbf{n}\).

\[\large{Autocorr(\mathbf{n}) = corr(series, \,serie\,con\,retardo\,\mathbf{n})}\]

Mide cuánto depende el valor de la serie en el presente del valor hace \(\mathbf{n}\) días.

El paquete feasts en R ofrece funciones para calcular y graficar la autocorrelación directamente desde objetos tsibble: - ACF() calculará la ACF. - autoplot() sobre el resultado de ACF() la graficará.

# Graficar ACF para data_monthly
gasoline_prod_monthly %>%
  ACF(value) %>%
  autoplot() +
  labs(title = "Función de Autocorrelación para Producción de Gasolina")

Los retardos están en el eje x. Ten en cuenta que la correlación con retraso-0 siempre es 1.

La banda azul punteada es el intervalo de confianza: si el valor de autocorrelación está fuera de ese rango, la dependencia es probablemente “real”, pero si está dentro, es probable que sea debido a la aleatoriedad.

Una advertencia: malinterpretación de la autocorrelación La autocorrelación parece ser muy fuerte: ¿eso significa que hay una dependencia muy fuerte entre los valores pasados y presentes que podemos usar para modelar? Echemos un vistazo a los datos nuevamente:

Estacionariedad

La mayoría de los modelos de series temporales asumen que la serie temporal es estacionaria. En lenguaje sencillo, esto significa que las características de la serie no dependen del tiempo en que ocurren (pero podrían depender de los valores pasados). Esto significa, entre otras cosas:

• Media constante
• Varianza constante

Las series temporales que presentan tendencia y estacionalidad no son estacionarias.

La mayoría de los modelos de series temporales asumen la estacionariedad de la serie que están modelando, por lo que a menudo debes eliminar la tendencia y la estacionalidad de la serie, construir el modelo y luego reintegrarlas. ## Eliminación de Tendencia

Existen algunas formas de eliminar la tendencia:

• Ajustar primero la curva de tendencia
• Diferenciación

En lugar de observar los valores reales de la serie, observa las diferencias en los valores:

\(U_t = Y_{t} - Y_{t-1}\)

La función difference() en R se puede utilizar para diferenciar la serie.

• Otras transformaciones

Los rendimientos y los logaritmos de los rendimientos, \(\log{\frac{Y_t - Y_{t-1}}{Y_t}}\), también se utilizan comúnmente en la modelización de precios de acciones.

---

Utilicemos la diferenciación para eliminar la tendencia de la serie de producción de gasolina.

# Diferenciar data_monthly
gasoline_prod_monthly_diff <- gasoline_prod_monthly %>%
  mutate(value_diff = difference(value))

head(gasoline_prod_monthly)

## # A tsibble: 6 x 2 [1M]
##   year_month value
##        <mth> <dbl>
## 1  1991 feb.  6.72
## 2  1991 mar.  7.09
## 3  1991 abr.  6.97
## 4  1991 may.  7.19
## 5  1991 jun.  7.62
## 6  1991 jul.  7.45

head(gasoline_prod_monthly_diff)

## # A tsibble: 6 x 3 [1M]
##   year_month value value_diff
##        <mth> <dbl>      <dbl>
## 1  1991 feb.  6.72     NA    
## 2  1991 mar.  7.09      0.371
## 3  1991 abr.  6.97     -0.113
## 4  1991 may.  7.19      0.218
## 5  1991 jun.  7.62      0.427
## 6  1991 jul.  7.45     -0.166

# Graficar los valores diferenciados
gasoline_prod_monthly_diff %>%
  autoplot(value_diff) +
  labs(
    title = "Producción de Gasolina Diferenciada Mensualmente",
    x = "Fecha",
    y = "Diferencia"
  )

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

# Graficar ACF para los datos diferenciados
gasoline_prod_monthly_diff %>%
  ACF(value_diff) %>%
  autoplot() +
  labs(title = "ACF de la Serie Diferenciada")

## Lidiando con la estacionalidad

Lidiar con la estacionalidad es mucho más complicado. Algunas cosas que puedes hacer por ti mismo:

• Modelar dentro de un período
• Usar datos ajustados estacionalmente
• Realizar el ajuste estacional tú mismo
  • Usar diferenciación de orden superior (diferenciación estacional)
  • Usar los residuos de STL() que están en el componente remainder.
• Usar modelos especiales como el Modelo ARIMA Estacional (SARIMA)

# Extraer residuos y eliminar valores NA
residuals_tsibble <- gasoline_prod_monthly %>%
  model(STL(value)) %>%
  components() %>%
  select(remainder) %>% 
  na.omit()

head(residuals_tsibble)

## # A tsibble: 6 x 2 [1M]
##   remainder year_month
##       <dbl>      <mth>
## 1  -0.0787   1991 feb.
## 2   0.0757   1991 mar.
## 3  -0.103    1991 abr.
## 4  -0.0210   1991 may.
## 5   0.175    1991 jun.
## 6  -0.00293  1991 jul.

# Graficar ACF de los residuos
residuals_tsibble %>%
  ACF(remainder) %>%
  autoplot() +
  labs(title = "ACF de los Residuos")

Tarea

1. Usa lag() y cor() para implementar la función de autocorrelación por ti mismo.

2. Es bastante fácil identificar cuándo una serie no es estacionaria si tiene tendencia o estacionalidad. Sin embargo, es más difícil afirmar que una serie que parece estacionaria realmente es estacionaria. Aprende sobre la prueba de Dickey-Fuller y la función adf.test() del paquete tseries.

La prueba de Dickey-Fuller aumentada (ADF) se utiliza para probar si una serie temporal es estacionaria o no.

Notas:

La hipótesis nula de la prueba de Dickey-Fuller aumentada es que hay una raíz unitaria (lo que implica no estacionariedad). La hipótesis alternativa es que no hay una raíz unitaria (estacionariedad). Si el valor p está por encima de un umbral crítico (generalmente 0.05), entonces no podemos rechazar la hipótesis nula de que hay una raíz unitaria.

# Realizar la prueba de Dickey-Fuller en la serie mensual
adf_result <- adf.test(gasoline_prod_monthly$value, alternative = "stationary")

# Imprimir los resultados de la prueba
print('Resultado de la prueba Dickey-Fuller:')

## [1] "Resultado de la prueba Dickey-Fuller:"

print(adf_result)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  gasoline_prod_monthly$value
## Dickey-Fuller = -2.4632, Lag order = 6, p-value = 0.3811
## alternative hypothesis: stationary

\[H_0: \text{La serie tiene una raíz unitaria (No Estacionaria)}\] \[VS\] \[ H_a:\text{la serie no tiene una raiz unitaria (Estacionaria)} \]

# Realizar la prueba de Dickey-Fuller en los datos diferenciados
adf_result_diff <- adf.test(na.omit(gasoline_prod_monthly_diff$value_diff), alternative = "stationary")

## Warning in adf.test(na.omit(gasoline_prod_monthly_diff$value_diff), alternative
## = "stationary"): p-value smaller than printed p-value

# Imprimir los resultados de la prueba
print('Resultado de la prueba Dickey-Fuller (Serie Diferenciada):')

## [1] "Resultado de la prueba Dickey-Fuller (Serie Diferenciada):"

print(adf_result_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(gasoline_prod_monthly_diff$value_diff)
## Dickey-Fuller = -12.63, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

En la práctica, buscamos que los datos sean estacionarios para poder proseguir con el tratamiento estadístico de los mismos.

La estacionariedad es un concepto fundamental para el análisis de las series de tiempo.

• La estacionariedad significa que el proceso mediante el cual se crearon los datos, es constante en el tiempo.
• CUIDADO: Esto no significa que los datos no puedan cambiar. Significa que los supuestos distribucionales del mecanismo que genera los datos son constantes en el tiempo.
• El punto de partida de la metodología Box-Jenkins (Construcción ARIMA) es la estacionariedad. De ahí su importancia práctica.

Modelos de series temporales

Dos modelos básicos

1. Ruido blanco

El mejor modelo para un valor futuro es simplemente una media constante.

\[\large{valor\_hoy = media + ruido\_hoy}\]

\[Y_t = \mu + \varepsilon_t\]

2. Camino aleatorio

El mejor modelo para un valor futuro es el valor anterior.

\[\large{valor\_hoy = valor\_ayer + ruido\_hoy}\]
\[Y_t = Y_{t-1} + \varepsilon_t\]

Modelo de promedio móvil (MA)

\[MA(q): \large{valor\_hoy = media + promedio(pasado \, \mathrm{q} \,ruidos)}\]

\[Y_t = \mu + \varepsilon_t + \phi_1 \varepsilon_{t-1} + \phi_2 \varepsilon_{t-2} +...+ \phi_q \varepsilon_{t-q} \]

Detección del orden del modelo MA

La ACF de un proceso \(MA(q)\) “se corta” después del rezago \(q\). Usaremos la ACF de la serie diferenciada para elegir el orden \(q\).

# ACF de los datos diferenciados (de nuevo para referencia)
gasoline_prod_monthly_diff %>%
  ACF(value_diff) %>%
  autoplot() +
  labs(title = "ACF de la Serie de Gasolina Diferenciada")

Aunque el patrón no es tan claro como en el ejemplo idealizado, podríamos decir que después del rezago 5, la autocorrelación es básicamente inexistente.

Por lo tanto, ajustaremos un modelo \(MA(5)\).

Ajuste del modelo MA(q)

Usaremos ARIMA() del paquete fable para especificar y ajustar el modelo.

1. Crear el modelo: modelo <- data %>% model(ARIMA(value ~ 0 + ma(q)))
2. Ajustar y resumir: report(modelo)

# Ajustar el modelo MA(5) a los datos mensuales
# Usamos d=1 para que el modelo trabaje con la serie diferenciada internamente
ma_model_fit <- gasoline_prod_monthly %>%
  model(ma5 = ARIMA(value ~ 0 + ma(1:5)))

## Warning: 1 error encountered for ma5
## [1] could not find function "ma"

# Imprimir resumen del modelo
report(ma_model_fit)

## Series: value 
## Model: NULL model 
## NULL model

Predicción con el modelo MA(q)

• forecast(modelo, h = ...) generará pronósticos para h períodos futuros.
• autoplot() se puede usar sobre el resultado para graficar las predicciones y los intervalos de confianza.

# Generar pronósticos desde el final de la serie
ma_forecast <- forecast(ma_model_fit, h = "4 years")

# Mostrar las predicciones
ma_forecast

## # A fable: 48 x 4 [1M]
## # Key:     .model [1]
##    .model year_month  value .mean
##    <chr>       <mth> <dist> <dbl>
##  1 ma5     2017 feb.     NA    NA
##  2 ma5     2017 mar.     NA    NA
##  3 ma5     2017 abr.     NA    NA
##  4 ma5     2017 may.     NA    NA
##  5 ma5     2017 jun.     NA    NA
##  6 ma5     2017 jul.     NA    NA
##  7 ma5     2017 ago.     NA    NA
##  8 ma5    2017 sept.     NA    NA
##  9 ma5     2017 oct.     NA    NA
## 10 ma5     2017 nov.     NA    NA
## # ℹ 38 more rows

# Graficar predicciones junto con los datos históricos
ma_forecast %>%
  autoplot(gasoline_prod_monthly %>% filter_index("2013-01-01" ~ .)) +
  labs(
    title = "Pronóstico de Producción de Gasolina con Modelo MA(5)",
    y = "Millones de Barriles"
  )

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning: Removed 48 rows containing missing values or values outside the scale range
## (`geom_line()`).

Nuevo conjunto de datos

En esta sección, trabajaremos con el conjunto de datos de desempleo de EE. UU. (US unemployment), que contiene datos trimestrales sobre el % de personas desempleadas en EE. UU. Los datos ya han sido ajustados estacionalmente.

# Cargar y preparar los datos de desempleo como tsibble
unemployment_url <- "https://raw.githubusercontent.com/datacamp/time-series-analysis-in-python-live-training/master/data/US_unemployment.csv"
unemployment_tsibble <- read_csv(unemployment_url) %>%
  mutate(date = yearquarter(as_date(date))) %>%
  as_tsibble(index = date)

## Rows: 203 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (1): value
## date (1): date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Graficar los datos
unemployment_tsibble %>%
  autoplot(value) +
  labs(
    title = "Tasa de Desempleo en EE. UU. (Ajustada Estacionalmente)",
    x = "Fecha", y = "% Desempleado"
  )

# Diferenciar los datos
unemployment_diff <- unemployment_tsibble %>%
  mutate(value_diff = difference(value))

# Graficar los datos diferenciados
unemployment_diff %>%
  autoplot(value_diff) +
  labs(
    title = "Tasa de Desempleo Diferenciada",
    x = "Fecha", y = "Diferencia"
  )

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

## Modelo autorregresivo (AR)

\[AR(p): \large{valor\_hoy = media + promedio(valores \, pasados \, \mathrm{p}) + ruido\_hoy}\]

\[Y_t = \mu + \theta_1 Y_{t-1} + \theta_2 Y_{t-2} + \dots + \theta_p Y_{t-p} + \varepsilon_t\]

Detección del orden del modelo

Para las series temporales que se comportan como un proceso \(AR(p)\), la Función de Autocorrelación Parcial (PACF) “se corta” en un cierto rezago. Ese rezago puede servir como el orden del modelo.

Aquí tienes la PACF de un verdadero modelo \(AR(2)\):

PACF de modelo AR(2)

# Graficar PACF para unemployment_diff
unemployment_diff %>%
  PACF(value_diff) %>%
  autoplot() +
  labs(title = "PACF para Desempleo Diferenciado")

La gráfica de la PACF sugiere que un modelo AR(2) podría ser apropiado, ya que hay picos significativos en los rezagos 1 y 2, y luego se corta.

Ajuste del modelo AR(p)

Usamos ARIMA() de fable con la especificación ar(p).

# Ajustar un modelo ARIMA(2,1,0) - un AR(2) sobre la serie diferenciada
ar_model_fit <- unemployment_tsibble %>%
  model(ar2_diff = ARIMA(value ~ 1 + ar(1:2), d = 1))

## Warning: 1 error encountered for ar2_diff
## [1] invalid type (list) for variable 'ar(1:2)'

# Imprimir resumen del modelo
report(ar_model_fit)

## Series: value 
## Model: NULL model 
## NULL model

\[\hat{y}_t' = c + \phi_1 y'_{t-1} + \phi_2 y'_{t-2}\] donde \(y'_t = y_t - y_{t-1}\)

\[\hat{desempleo}'_t = 0.0075 + 0.2552 \times desempleo'_{t-1} + 0.1702 \times desempleo'_{t-2}\]

# Pronóstico con el modelo AR(2)
ar_forecast <- forecast(ar_model_fit, h = "7 years") # 28 trimestres

# Graficar pronósticos
ar_forecast %>%
  autoplot(unemployment_tsibble) +
  labs(
    title = "Pronóstico de Tasa de Desempleo con Modelo ARIMA(2,1,0)",
    y = "% Desempleo"
  )

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning: Removed 28 rows containing missing values or values outside the scale range
## (`geom_line()`).

Proceso	Función de autocorrelación (ACF)	Función de autocorrelación parcial (ACFP)
MA(q)	Solo los q primeros coeficientes son significativos. El resto se anulan bruscamente (coef. 0 para retardo >q)	Decrecimiento rápido exponencial atenuado u ondas sinusoidales.
AR(p)	Decrecimiento rápido exponencial atenuado u ondas sinusoidales.	Solo los p primeros coeficientes son significativos. El resto se anulan bruscamente (coef. 0 para retardo >q)
ARIMA(p, d, q)	Comportamiento irregular en los retardos(1, … , q) con q picos. Decrecimiento para retardos posteriores a q.	Decrece (aproximadamente con exponenciales atenuados y ondas sinusoidales). No cero pronto.