En la Clase 3 pronosticamos series de tiempo utilizando métodos ingenuos y regresiones lineales. Si bien el modelo de regresión lineal no fue diseñado para manipular y pronosticar series de tiempo, es posible generar características1 y transformar el problema de pronóstico en un problema de regresión. Sin embargo, el supuesto para ajustar un modelo de regresión por mínimos cuadrados ordinarios (MCO)2 acerca del término del error \(\epsilon\) (independiente e idénticamente distribuido con media \(0\) y varianza \(\sigma^2\)) en muchos casos es violado al trabajar con series de tiempo. Existe la posibilidad de que los residuales del modelo aún presenten correlación3, por lo que necesitamos disponer de modelos que aborden dichas dependencias en el tiempo de manera natural.
En esta clase estudiaremos modelos univariados para pronóstico de series de tiempo, revisando los principios teóricos y prácticos de cada uno de ellos utilizando ecosistemas de R. Luego de esta clase podrás:
- Utilizar modelos de suavizamiento exponencial (Exponential Smoothing Models) para pronosticar series de tiempo
- Utilizar modelos autoregresivos (
AR), de medias móviles (MA), integrados (ARIMA) y estacionales (SARIMA) para pronosticar series de tiempo. - Modelar residuales de una regresión lineal con modelos de series de tiempo
- Abordar estacionalidades complejas con modelos adecuados (
batsytbats) - Utilizar modelos de regresión no-lineal como Prophet
Suavizamiento Exponencial (Exponential Smoothing)4 📈
El Suavizamiento Exponencial fue propuesto a finales de los 50’s (Brown, Holt, Winters) y ha motivado el desarrollo de métodos de pronósticos exitosos en el dominio de series de tiempo. Los pronósticos producidos por métodos de suavizamiento exponencial son medias ponderadas de observaciones pasadas, con pesos5 que decaen de manera exponencial a medida que las observaciones se hacen más distantes en el tiempo. Es decir, las observaciones recientes poseen mayor relevancia para efectos de elaboración de una predicción para valores futuros.
Suavizamiento Exponencial Simple (SES)
El Suavizamiento Exponencial Simple es un método adecuado para pronosticar datos que no poseen un patrón de tendencia o estacionalidad claro, siendo un método simple y rápido de computar. SES supone que la serie de tiempo sólo contiene un componente de un nivel \(\ell_t\) (level), que se mantiene en el tiempo y un error. Podemos representar el suavizamiento exponencial simple utilizando ecuaciones de pronóstico y suavizamiento de la siguiente manera6:
\[ \begin{align*} \text{Ecuación de Pronóstico} && \hat{y}_{t+h|t} & = \ell_{t}\\ \text{Ecuación de Suavizamiento} && \ell_{t} & = \alpha y_{t} + (1 - \alpha)\ell_{t-1}, \end{align*} \]
Podemos ver que el pronóstico para SES corresponde al nivel estimado para la observación más reciente en período \(t\). Para obtener el nivel \(\ell_t\), utilizamos la ecuación de suavizamiento, que actualiza el nivel anterior \(\ell_{t-1}\) mediante la integración de la información más reciente asociada a la observación \(y_t\). En este caso, \(\alpha\) es el ponderador o la constante de suavizamiento que toma valores entre \(0 \leq \alpha \leq 1\). El algoritmo aprende el nuevo nivel a partir de los datos observados más recientes. Debemos notar que es necesario inicializar el algoritmo definiendo una condición inicial para \(\ell_1\)7. Si sustituimos \(\ell_t\) en su fórmula obtenemos:
\[\begin{align*} \ell_{t} & = \alpha y_t + (1-\alpha) \left[\alpha y_{t-1} + (1-\alpha) \ell_{t-2}\right] \\ & = \alpha y_t + \alpha(1-\alpha) y_{t-1} + (1-\alpha)^2 \ell_{t-2} \\ & = \dots \\ & = \alpha y_t + \alpha(1-\alpha) y_{t-1} + \alpha(1-\alpha)^2 y_{t-2} + \dots \end{align*}\]
en este caso, terminamos con un promedio ponderado de todos los valores de la serie con pesos que decaen de manera exponencial8.
La constante de suavizamiento \(\alpha\) determina qué tanto peso9 se otorga a una observación pasada:
- \(\alpha = 1\): observaciones pasadas no influyen en los pronósticos10
- \(\alpha = 0\): observaciones pasadas poseen una influencia equitativa en los pronósticos
En la literatura, el parámetro \(\alpha\) posee como valores típicos 0.1 o 0.211 . Sin embargo, es posible obtener su valor óptimo minimizando alguna métrica de desempeño de interés, con el cuidado de no sobre-ajustar el método a los datos.
Finalmente, es posible entender SES como un proceso de aprendizaje adaptativo si consideramos la ecuación del nivel de la siguiente forma:
\[ \ell_{t} = \alpha y_t + (1-\alpha)\ell_{t-1} = \alpha y_t + \ell_{t-1} - \alpha \ell_{t-1} \]
así, podemos sustituir \(\ell_{t}\) por \(\hat{y}_{t+h|t}\)
\[\begin{align} \hat{y}_{t+1|t} & = \ell_{t} = \ell_{t-1} + \alpha(y_{t} - \ell_{t-1}) \\ & = \hat{y}_{t+1|t} + \alpha(y_t - \hat{y}_{t+1|t}) \\ & = \hat{y}_{t+1|t} + \alpha \varepsilon_t \end{align}\]
y notamos que el pronóstico es una actualización del pronóstico anterior más un componente que depende del error de dicho pronóstico anterior.
Show code
# SES
model_spec_ses <- exp_smoothing( # SES
error = "additive",
trend = "none",
seasonal = "none"
) %>%
set_engine("ets")
model_fit_ses <- model_spec_ses %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_ses
parsnip model object
Fit time: 20ms
ETS(A,N,N)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.9999
Initial states:
l = 2025924.8735
sigma: 285805.8
AIC AICc BIC
5869.033 5869.160 5878.821 Show code
calibration_tbl <- model_fit_ses %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Métodos con Tendencia
Tendencia Lineal de Holt
En 1957, Holt extendió SES para permitir que los pronósticos pudiesen incorporar efectos de tendencia. Este método involucra una ecuación de pronóstico, una ecuación de nivel y otra de tenencia según:
\[\begin{align*} \text{Ecuación de pronóstico}&& \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} \\ \text{Ecuación de nivel} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ \text{Ecuación de tendencia} && b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1}, \end{align*}\]
donde \(b_t\) denota al estimador de la tendencia (pendiente) de la serie en el período \(t\) cuyo parámetro de suavizamiento es \(\beta^*\), que toma valores entre \(0 \leq \beta^* \leq 1\).
Tal como con SES, la ecuación de nivel muestra que \(\ell_t\) es un promedio ponderado de la observación \(y_t\) y el pronóstico a un paso para el período \(t\), dado por \(\ell_t + b_{t-1}\). La ecuación de la tendencia muestra que \(b_t\) es un promedio ponderado de la tendencia estimada para el período \(t\) basada en \(\ell_t - \ell_{t-1}\) y \(b_{t-1}\), el estimador previo de la tendencia. Es decir, la tendencia previa se actualiza utilizando la diferencia entre los valores del nivel más actuales. En ese sentido, podemos hablar de una tendencia local que se adapta a través del tiempo mediante el parámetro de velocidad \(\beta^*\).
En este modelo, el pronóstico ya no es plano, sino posee tendencia, donde el pronóstico a \(h\) pasos es equivalente al último nivel estimado más \(h\) veces el último valor de la tendencia estimada12.
Show code
# HOLT
model_spec_holt <- exp_smoothing( # SES
error = "additive",
trend = "additive",
seasonal = "none"
) %>%
set_engine("ets")
model_fit_holt <- model_spec_holt %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_holt
parsnip model object
Fit time: 20ms
ETS(A,A,N)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.9999
beta = 0.0029
Initial states:
l = 2460121.725
b = -24350.0269
sigma: 284956.9
AIC AICc BIC
5869.853 5870.174 5886.167 Show code
calibration_tbl <- model_fit_holt %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Tendencia amortiguada
Como pudimos ver, los pronósticos generados por el método de Holt muestran una tendencia constante (creciente o decreciente) hacia el futuro, lo que a largo plazo puede tender a sobreestimar valores futuros. La introducción de un parámetro \(\phi\) para amortiguar (damping)13 la tendencia a una línea plana en el futuro:
\[\begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t} \\ \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)\phi b_{t-1}. \end{align*}\]
donde \(0 < \phi < 1\). Si \(\phi = 1\) obtenemos el método lineal de Holt.
Show code
# HOLT DAMPED
model_spec_holt_damped <- exp_smoothing( # SES
error = "additive",
trend = "additive",
seasonal = "none",
damping = "damped"
) %>%
set_engine("ets")
model_fit_holt_damped <- model_spec_holt_damped %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_holt_damped
parsnip model object
Fit time: 15ms
ETS(A,Ad,N)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.9999
beta = 1e-04
phi = 0.8236
Initial states:
l = 2460122.2229
b = -44350.7438
sigma: 284326.7
AIC AICc BIC
5869.975 5870.426 5889.551 Show code
calibration_tbl <- model_fit_holt_damped %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Lamentablemente, para la serie de USgas el modelo no presenta mayor aporte, sin embargo, podemos ver un ejemplo de aplicación para pronosticar la población australiana acá.
Métodos con Estacionalidad
Consideremos los tipos de tendencia y estacionalidad del siguiente video. En el caso de tendencia aditiva, la tendencia crece/decrece de manera lineal a través del tiempo. La tendencia multiplicativa implica crecimiento/decrecimiento exponencial en el tiempo. Finalmente, la tendencia polinomial permite capturar una variedad de formas de manera flexible. De manera similar, hallamos estacionalidad aditiva, donde cada estación difere en una cantidad específica de otra y estacionalidad multiplicativa, donde cada estación difere de otra en un porcentaje específico.
Podemos entender un modelo aditivo como la suma de los componentes sistemáticos de una serie más el componente de ruido o no-sistemático de la siguiente forma:
\[y_t = \text{nivel}\ + \text{tendencia}\ + \text{estacionalidad}\ + \text{ruido}\]
y a los modelos multiplicativos según:
\[y_t = \text{nivel}\ \times \text{tendencia}\ \times \text{estacionalidad}\ \times \text{ruido}\]
En 1960 Holt y Winters extendieron el método de Holt para capturar efectos estacionales. El método estacional de Holt-Winters involucra una ecuación de pronóstico y tres ecuaciones de suavizamiento para el nivel \(\ell_t\), la tendencia \(b_t\) y la estacionalidad \(s_t\), con parámetros \(\alpha\), \(\beta*\) y \(\gamma\), respectivamente. En este caso, usamos \(m\) para denotar el período asociado a la estacionalidad.
Holt-Winters Aditivo
Cuando las variaciones estacionales son relativamente constantes a través de la serie es preferible utilizar el modelo de Holt-Winters Aditivo. En este caso, la componente estacional se expresa en términos absolutos según la escala de la serie observada, y la ecuación del nivel es ajustada estacionalmente mediante la sustracción del componente estacional, como sigue:
\[\begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} + s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha(y_{t} - s_{t-m}) + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1}\\ s_{t} &= \gamma (y_{t}-\ell_{t-1}-b_{t-1}) + (1-\gamma)s_{t-m}, \end{align*}\]
donde \(k\) es la parte entera de \((h-1)/m\)14. La ecuación del nivel presenta un promedio ponderado entre la observación ajustada por estacionalidad y el pronóstico no-estacional (\(\ell_{t-1}+b_{t-1}\)) para el período \(t\). La ecuación de la tendencia es la misma que para el método de Holt. LA ecuación de estacionalidad muestra un promedio ponderado entre el indice estacional actual, , y el índice estacional de la misma temporada en el último período15.
Show code
# HOLT WINTERS aditivo
model_spec_hw <- exp_smoothing(
error = "additive",
trend = "additive",
season = "additive",
seasonal_period = 12
) %>%
set_engine("ets")
model_fit_hw <- model_spec_hw %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_hw
parsnip model object
Fit time: 131ms
ETS(A,A,A)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.4874
beta = 1e-04
gamma = 1e-04
Initial states:
l = 1878651.3185
b = 3462.6919
s = 458272 -47106.28 -289819.3 -385374.7 -181214.9 -199612.4
-358117.1 -344574.9 -157620.9 269526 491002.5 744640
sigma: 100392.1
AIC AICc BIC
5478.497 5481.994 5533.962 Show code
calibration_tbl <- model_fit_hw %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Holt-Winters Multiplicativo
Por otra parte, el modelo multiplicativo es preferido cuando las variaciones estacionales cambian de manera proporcional con el nivel de la serie. En este caso, la componente estacional se expresa en términos relativos (porcentajes), y la serie es ajustada estacionalmente mediante la división a través del componente estacional.
\[\begin{align*} \hat{y}_{t+h|t} &= (\ell_{t} + hb_{t})s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha \frac{y_{t}}{s_{t-m}} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ b_{t} &= \beta^*(\ell_{t}-\ell_{t-1}) + (1 - \beta^*)b_{t-1} \\ s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + b_{t-1})} + (1 - \gamma)s_{t-m}. \end{align*}\]
Show code
# HOLT WINTERS multiplicativo
model_spec_hw_m <- exp_smoothing(
error = "multiplicative",
trend = "multiplicative",
season = "multiplicative",
seasonal_period = 12
) %>%
set_engine("ets")
model_fit_hw_m <- model_spec_hw_m %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_hw_m
parsnip model object
Fit time: 235ms
ETS(M,Md,M)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.5415
beta = 1e-04
gamma = 1e-04
phi = 0.9768
Initial states:
l = 1848879.3258
b = 1.0017
s = 1.2334 0.9702 0.8501 0.8058 0.9087 0.8969
0.8191 0.8289 0.9239 1.1388 1.2468 1.3774
sigma: 0.0467
AIC AICc BIC
5444.938 5448.869 5503.666 Show code
calibration_tbl <- model_fit_hw_m %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Holt-Winters Amortiguado
Es posible amortigüar los métodos de Holt-Winters (aditivos y multiplicativos) revisados anteriormente. Un ejemplo es el método de Holt-Winters con tendencia amortiguada y estacionalidad multiplicativa que sigue:
\[\begin{align*} \hat{y}_{t+h|t} &= \left[\ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t}\right]s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha(y_{t} / s_{t-m}) + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)\phi b_{t-1} \\ s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + \phi b_{t-1})} + (1 - \gamma)s_{t-m}. \end{align*}\]
Show code
# HOLT WINTERS multiplicativo con damping
model_spec_hw_m_damped <- exp_smoothing(
error = "multiplicative",
trend = "multiplicative",
season = "multiplicative",
damping = "damped",
seasonal_period = 12
) %>%
set_engine("ets")
model_fit_hw_m_damped <- model_spec_hw_m_damped %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_hw_m_damped
parsnip model object
Fit time: 136ms
ETS(M,Md,M)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets,
Call:
alpha = alpha, beta = beta, gamma = gamma)
Smoothing parameters:
alpha = 0.5415
beta = 1e-04
gamma = 1e-04
phi = 0.9768
Initial states:
l = 1848879.3258
b = 1.0017
s = 1.2334 0.9702 0.8501 0.8058 0.9087 0.8969
0.8191 0.8289 0.9239 1.1388 1.2468 1.3774
sigma: 0.0467
AIC AICc BIC
5444.938 5448.869 5503.666 Show code
calibration_tbl <- model_fit_hw_m_damped %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
Modelos de espacio de estado (ETS)
Las combinaciones y variaciones de componentes de tendencia y estacionalidad hacen posible nueve métodos de suavizamiento exponencial, denotados por el par \((T,S)\) según el tipo de tendencia y estacionalidad (Trend, Seasonality). Así el método \((A_d,M)\) es un suavizamiento exponencial con tendencia amortiguada y estacionalidad multiplicativa.
Show code
knitr::include_graphics("https://otexts.com/fpp3/figs/pegelstable-1.png")
](https://otexts.com/fpp3/figs/pegelstable-1.png)
Figure 1: Fuente: A taxonomy of exponential smoothing methods (fpp3)
Los métodos de suavizamiento exponencial poseen modelos estadísticos para generar intervalos de predicción16 y producir una distribución de pronóstico. Cada modelo consiste de una ecuación de medición (measurement equation) que describe los datos observados y algunas ecuaciones de estado (state equations) que describen la manera en que los componentes no observados o estados (nivel, tendencia, estacionalidad) cambian en el tiempo. De ahí que se conozcan como modelos de espacio de estados (state space models).
Para distinguir entre un modelo con errores aditivos de uno con errores multiplicativos, se agrega una tercera letra a la clasificación anterior, etiquetando a cada modelo de espacio de estados como ETS(Error, Trend, Seasonal). Podemos ver todos los modelos de espacio de estado en el siguiente cuadro:
Show code
knitr::include_graphics("https://otexts.com/fpp3/figs/statespacemodels-1.png")
](https://otexts.com/fpp3/figs/statespacemodels-1.png)
Figure 2: Fuente: Innovations state space models for exponential smoothing (fpp3)
Consideremos el suavizamiento exponencial simple con errores aditivos: ETS(A,N,N)
\[\begin{align*} \text{Ecuación de pronóstico} && \hat{y}_{t+1|t} & = \ell_{t}\\ \text{Ecuación de suavizamiento} && \ell_{t} & = \alpha y_{t} + (1 - \alpha)\ell_{t-1}. \end{align*}\]
Si se reordena la ecuación de suavizamiento para el nivel, se obtiene la forma de corrección de error que sigue:
\[\begin{align*} \ell_{t} %&= \alpha y_{t}+\ell_{t-1}-\alpha\ell_{t-1}\\ &= \ell_{t-1}+\alpha( y_{t}-\ell_{t-1})\\ &= \ell_{t-1}+\alpha e_{t}, \end{align*}\]
donde \(e_{t}=y_{t}-\ell_{t-1}=y_{t}-\hat{y}_{t|t-1}\) es el residual del período \(t\). Por lo tanto, podemos escribir \(y_t = \ell_{t-1}+ e_{t}\), tal que cada observación pueda ser representada por un nivel previo más un error. Si se asume que los residuales (errores de entrenamiento a un paso) \(e_t\) distribuyen como ruido blanco gaussiano, entonces las ecuaciones pueden escribirse:
\[\begin{align} y_t &= \ell_{t-1} + \varepsilon_t\\ \ell_t&=\ell_{t-1}+\alpha \varepsilon_t \end{align}\]
Donde la primera ecuación corresponde a la ecuación de medición u observación y la segunda corresponde a la ecuación de estado o transición. Estas ecuaciones, junto a la distribución estadística del error, forman un modelo estadístico conocido como modelo de espacio de estado de innovaciones (innovations state space model) subyacente al modelo de suavizamiento exponencial simple asociado. Todas las ecuaciones utilizan las innovaciones provienen del mismo proceso de error aleatorio, \(\varepsilon_t\), por lo que también son conocidos como modelos de fuente única de error.
Modelos ARIMA 📈
Remember that all models are wrong; the practical question is how wrong do they have to be to not be useful1718
George E. P. Box
En la sección anterior discutimos técnicas de pronóstico que, en términos generales, se basan en la idea de representar la serie de tiempo como una suma de componentes distintos: uno determinístico y otro estocástico19. El componente determinístico se modela mediante una función del tiempo, mientras que para el componente estocástico se asume un ruido blanco que es adicionado a la señal determinística generando un comportamiento estocástico en la serie de tiempo. Un supuesto importante es que el ruido aleatorio es generado por un proceso de elaboración de innovaciones o shocks independientes. Como se ha mencionado anteriormente, en la práctica dicho supuesto usualmente es violado, por lo que las observaciones sucesivas presentan dependencias seriales para las cuales un modelo de suavizamiento exponencial pueda ser ineficiente e incluso inapropiado20. Para incorporar de manera formal esta estructura de dependencia, a continuación exploraremos una clase general de modelos conocidos como modelos autoregresivos integrados de medias móviles (Autoregressive Integrated Moving Average) o ARIMA21.
Diferenciación
En la Clase 2 introdujimos el concepto de estacionariedad para describir una serie de tiempo cuyas propiedades estadísticas no dependen del momento en que la serie es observada. Por lo tanto, series con tendencia o estacionalidad no son estacionarias, pues la tendencia y la estacionalidad afectarían el valor de la serie de tiempo en diferentes momentos22. Esto último puede ser confuso: una serie de tiempo puede exhibir un comportamiento cíclico (aunque sin tendencia ni estacionalidad) y ser estacionaria, ya que los ciclos no poseen un largo fijo. En general una serie estacionaria no presentará patrones predecibles al largo plazo. El gráfico debería parecer una serie relativamente horizontal, con algunas variaciones cíclicas y varianza constante.
Revisemos el precio de las acciones de Google:
Show code
quantmod::getSymbols("GOOGL")
[1] "GOOGL"En la figura anterior podemos ver el precio de cierre de las acciones de Google, una serie claramente no estacionaria. Sin embargo, si calculamos la diferencia \((y_t - y_{t-1})\) entre sus observaciones consecutivas podemos ver que la serie ahora es aparentemente estacionaria:
Transformaciones tales como logarítmos (o Box-Cox) pueden ayudar a estabilizar la varianza de una serie de tiempo. En este caso, la diferenciación logra estabilizar la media de una serie de tiempo removiendo los cambios en el nivel de una serie, y por ende, eliminado (o reduciendo) su tendencia y estacionalidad. Tal como vimos en clases previas, el gráfico de la ACF es útil para identificar series de tiempo no estacionarias23.
La ACF de la diferencia de precios asociada a las acciones de Apple se ve idéntica la del ruido blanco. El test de Ljung-Box (p-value = 0.6257) sugiere que el cambio diario del precio de las acciones de Google es esencialmente un valor aleatorio que no está correlacionado con días previos.
Show code
Box-Ljung test
data: .
X-squared = 8.0323, df = 10, p-value = 0.6257Caminata Aleatoria (Random Walk)
La serie diferenciada es el cambio entre observaciones consecutivas de la serie original y puede ser escrito como:
\[y'_t = y_t - y_{t-1}\]
La serie diferenciada tendrá \(T-1\) valores. Si dichos valores son ruido blanco, el modelo para la serie original puede ser escrito como
\[y_t - y_{t-1} = \varepsilon_t\]
donde \(\varepsilon_t\) denota el ruido blanco. Reordenando, obtenemos el modelo de caminata aleatoria (random walk) 😓
\[y_t = y_{t-1} + \varepsilon_t\]
Las caminatas aleatorias son ampliamente utilizadas para datos no-estacionarios, particularmente financieros o económicos, ya que típicamente poseen:
- Períodos largos de tendencia al crecimiento/decrecimiento aparente
- Cambios repentinos y poco predecibles en su dirección
Los pronósticos asociados a las caminatas aleatorias son equivalentes a la última observación24, dado que los movimientos futuros son impredecibles, con la misma probabilidad de ascender o descender.
En particular, si la diferencia no posee media cero, entonces
\[y_t - y_{t-1} = c + \varepsilon_t\quad\text{o}\quad {y_t = c + y_{t-1} + \varepsilon_t}\]
lo que se conoce como caminata aleatoria con deriva (random walk with drift).
Show code
# Shumway & Stoffer example
set.seed(154)
w <- rnorm(200)
x <- cumsum(w)
wd <- w +.2 # drift
xd <- cumsum(wd)
rw_plot <- tibble(
index = 1:length(w),
x,
xd
) %>%
pivot_longer(
cols = -index,
names_to = "series"
) %>%
ggplot(aes(x = index, y = value, color = series)) +
geom_line() +
geom_abline(intercept = 0, slope = 0, linetype="dotted") +
geom_abline(intercept = 0, slope = 0.2, linetype="dotted") +
labs(title = "Random Walk y Random Walk w/drift") +
theme_bw()
plotly::ggplotly(rw_plot)
Diferenciación estacional
El mismo principio de diferenciación aplica en series estacionales, donde la diferencia entra la observación y la observación previa de la misma temporada se conoce como diferencia estacional25. En ocasiones será necesario aplicar diferencias estacionales y primeras diferencias para obtener datos estacionarios, existiendo un grado de subjetividad al momento de aplicar diferencias suficientes para obtener datos estacionarios. Debido a lo anterior, es importante que las diferencias sean interpretables. Para determinar de manera más objetiva si una serie requiere de diferencias o no es posible realizar el test de raíz unitaria26.
Podemos ilustrar el proceso iterativo de diferenciación con la serie estacional AirPassengers
y sus respectivas funciones de autocorrelación
Show code
AirPassengers %>%
timetk::tk_tbl() %>%
timetk::plot_acf_diagnostics(
.date_var = index,
.value = value,
.lags = 24,
.interactive = TRUE
)
Modelos Autoregresivos (AR)
Consideremos un modelo recursivo en el cual la observación de hoy depende de la de ayer para todo momento \(t\), hoy = constante + pendiente x ayer + ruido. Dicho proceso podría estar centrado, cuya relación sería (hoy - media) = pendiente x (ayer - media) + ruido. De manera más formal, si \(\varepsilon_{t}\) es ruido blanco con media cero y varianza \(\sigma^2\), entonces
\[y_{t} - \mu = \phi(y_{t-1} - \mu) + \varepsilon_t\]
se considera un proceso autoregresivo de orden 1, AR(1), si \(\phi \neq 0\). Valores grandes de \(\phi\) indicarán una alta autocorrelación, mientras que si \(\phi < 0\) la serie resultante será oscilatoria.
Simulemos tres procesos AR(1):
Show code
set.seed(2022)
sim_ar <- tibble::tibble(Index = 1:100,
ar_05 = arima.sim(model = list(ar=0.5), n = 100),
ar_09 = arima.sim(model = list(ar=0.9), n = 100),
ar_m075 = arima.sim(model = list(ar=-0.75), n = 100))
plot_sim_ar <- sim_ar %>%
tidyr::pivot_longer(cols = -(Index),
names_to = "Model",
values_to = "Series") %>%
ggplot(aes(x = Index, y = Series, color = Model)) +
geom_line() +
theme_bw()
plotly::ggplotly(plot_sim_ar)
cuyas respectivas funciones de autocorrelación son:
De forma más general, un modelo autoregresivo de orden \(p\), \(AR(p)\), tiene la forma
\[y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_{t}\]
e involucra la regresión de una variable consigo misma. Aquí \(y_t\) es estacionaria y \(\varepsilon_{t}\) es ruido blanco. Los parámetros \(\phi_{1}, \dots, \phi_{p}\) son constantes y la media de \(y_t\) es cero27. Note que para el modelo AR(1):
- Si \(\phi_1 = 0\) = 0 y \(c = 0\), \(y_t\) es ruido blanco
- Si \(\phi_1 = 1\) = 0 y \(c = 0\), \(y_t\) es una caminata aleatoria
- Si \(\phi_1 = 1\) = 0 y \(c \neq 0\), \(y_t\) es una caminata aleatoria con deriva
- Si \(\phi_1 < 0\), \(y_t\) oscila alrededor de la media.
Usualmente los parámetros \(\phi\) y sus combinaciones lineales son restringidos a oscilar entre [-1,1].
Modelos de Medias Móviles (MA)
Los modelos con autocorrelacion pueden ser construidos a partir del ruido blanco, un modelo que pondera el ruido pasado y actual es llamado media móvil, como por ejemplo hoy = media + ruido + pendiente x (ruido de ayer), lo que puede escribirse de manera más formal de la siguiente forma
\[y_{t} = \mu + \varepsilon_t + \theta\varepsilon_{t-1} \]
que corresponde a una media móvil de orden 1, MA(1). Simulemos tres procesos MA(1):
Show code
set.seed(2022)
sim_ma <- tibble::tibble(Index = 1:100,
ma_05 = arima.sim(model = list(ma=0.5), n = 100),
ma_09 = arima.sim(model = list(ma=0.9), n = 100),
ma_m05 = arima.sim(model = list(ma=-0.5), n = 100))
plot_sim_ma <- sim_ma %>%
tidyr::pivot_longer(cols = -(Index),
names_to = "Model",
values_to = "Series") %>%
ggplot(aes(x = Index, y = Series, color = Model)) +
geom_line() +
theme_bw()
plotly::ggplotly(plot_sim_ma)
cuyas respectivas funciones de autocorrelación son:
En el caso de las medias móviles, en lugar de utilizar valores pasados de la variable a pronosticar, se utilizan errores de pronóstico para realizar un modelo del tipo regresión. Nos referimos al modelo de medias móviles de orden \(q\), \(MA(q)\) con la expresión
\[y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q}\]
donde \(\varepsilon_{t}\) es ruido blanco. En este caso, cada valor de \(y_t\) puede pensarse como una media móvil ponderada de los errores de pronóstico pasados.
Modelos Autoregresivos de Medias Móviles (ARMA)
La unión de las aproximaciones anteriores en un sólo modelo permite generar una regresión múltiple con observaciones rezagadas y errores como predictores del siguiente modo
\[y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q} + \varepsilon_t\] que se conoce como el modelo ARMA(p,q) con las últimas p observaciones y últimos q errores como regresores. Los modelos ARMA sólo funcionan con datos estacionarios, por lo que si el proceso no es estacionario será necesario diferenciar previamente la serie. Si una serie necesita ser diferenciada d veces para hacerla estacionaria, el modelo resultante es llamado un Modelo Autoregresivo Integrado de Medias Móviles o ARIMA(p,d,q)
Modelos Autoregresivos Integrados de Medias Móviles (ARIMA)
Si se combinan los procesos de diferenciación con los modelos autoregresivos y de medias móviles, obtenemos el modelo autoregresivo integrado de medias móviles (ARIMA)28. En este contexto, con integración nos referimos al proceso inverso de la diferenciación.
\[\begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t} \end{equation}\]
donde \(y'_{t}\) es la serie diferenciada29. El modelo \(ARIMA(p,d,q)\) describe un modelo autoregresivo de orden \(p\), con medias móviles de orden \(q\) con grado de integración \(d\).
Si definimos el operador de diferencias \(B y_{t} = y_{t - 1}\), entonces \(y'_{t} = y_{t} - y_{t-1} = y_t - By_{t} = (1 - B)y_{t}\) y \(B(By_{t}) = B^{2}y_{t} = y_{t-2}\). De manera más general, la diferencia de orden \(d\) puede escribirse como \((1 - B)^{d} y_{t}\), luego:
\[\begin{equation} \begin{array}{c c c c} (1-\phi_1B - \cdots - \phi_p B^p) & (1-B)^d y_{t} &= &c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t\\ {\uparrow} & {\uparrow} & &{\uparrow}\\ \text{AR($p$)} & \text{$d$ diferencias} & & \text{MA($q$)}\\ \end{array} \end{equation}\]
Aplicando el modelo ARIMA a la serie us_monthly
Show code
parsnip model object
Fit time: 14.1s
Series: outcome
ARIMA(1,0,2)(0,1,1)[12] with drift
Coefficients:
ar1 ma1 ma2 sma1 drift
0.9799 -0.5114 -0.3233 -0.8763 2232.3015
s.e. 0.0263 0.0732 0.0686 0.0807 811.7828
sigma^2 estimated as 9.76e+09: log likelihood=-2344.13
AIC=4700.27 AICc=4700.75 BIC=4719.46Show code
calibration_tbl <- model_fit_arima %>%
# 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 %>% rmarkdown::paged_table()
Show code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
La implementación de auto.arima()30 es una algoritmo que funciona del siguiente modo:
- Selecciona la cantidad de diferencias
dmediante tests de raíz unitaria. - Selecciona
pyqminimizando \(AIC_{c}\). - Estima los parámetros del modelo mediante la estimación de su máxima verosimilitud.
- Busca a través del espacio de modelos por pasos (stepwise).
Modelos Autoregresivos Integrados de Medias Móviles Estacionales (SARIMA)
Revisemos el modelo ajustado anteriormente para us_monthly usando auto.arima()
Show code
model_fit_arima
parsnip model object
Fit time: 14.1s
Series: outcome
ARIMA(1,0,2)(0,1,1)[12] with drift
Coefficients:
ar1 ma1 ma2 sma1 drift
0.9799 -0.5114 -0.3233 -0.8763 2232.3015
s.e. 0.0263 0.0732 0.0686 0.0807 811.7828
sigma^2 estimated as 9.76e+09: log likelihood=-2344.13
AIC=4700.27 AICc=4700.75 BIC=4719.46el modelo ajustado posee estacionalidad en una media móvil, sma1. El modelo ajustado corresponde a un modelo ARIMA(1,0,2)(0,1,1)[12], es decir, un modelo ARIMA con componente estacional dada por la forma ARIMA(p,d,q)(P,D,Q)[m], donde:
des la cantidad de diferencias con respecto al rezago 1pes la cantidad de rezagos AR ordinarios: \(y_{t-1},y_{t-2},\dots,y_{t-p}\)qes la cantidad de rezagos MA ordinarios: \(\varepsilon_{t-1},\varepsilon_{t-2},\dots,\varepsilon_{t-q}\)Des la cantidad de diferencias estacionalesPes la cantidad de rezagos AR estacionales: \(y_{t-m},y_{t-2m},\dots,y_{t-Pm}\)Qes la cantidad de rezagos MA estacionales: \(\varepsilon_{t-m},\varepsilon_{t-2m},\dots,\varepsilon_{t-Qm}\)mes la cantidad de observaciones por ciclo31
Modelos de Regresión Dinámica 📈
Hasta ahora hemos revisado modelos de pronóstico que utilizan sólo la historia de la serie, sin considerar información adicional que podría estar disponible. Pero a menudo puede existir información adicional32 que nos permita mejorar las predicciones.
Los Modelos de Regresión Dinámica (Dynamic Regression Models) corresponden a un modelo de regresión lineal con errores ARIMA:
\[y_t = \beta_0 + \beta_{1}x_{1,t} + \dots + \beta_{r}x_{r,t} + \varepsilon_t\]
donde \(y_t\) es modelado como función de \(r\) variables independientes y el término del error es un proceso ARIMA (en lugar del ruido blanco usual).
Show code
# Regresion Lineal con errores ARIMA
model_fit_linreg_arima <- model_spec_arima %>%
fit(
y ~ date
+ lubridate::month(date, label = TRUE)
+ fourier_vec(date, period = 12)
+ fourier_vec(date, period = 24)
+ fourier_vec(date, period = 48),
data = training(splits)
)
model_fit_linreg_arima
parsnip model object
Fit time: 6.6s
Series: outcome
Regression with ARIMA(1,1,1)(0,0,2)[12] errors
Coefficients:
ar1 ma1 sma1 sma2 drift
0.4036 -0.9150 0.0651 -0.2282 2237.8383
s.e. 0.0762 0.0301 0.0777 0.0810 828.4225
fourier_vec_date_period_12 fourier_vec_date_period_24
-50147.73 1487.232
s.e. 23403.87 10829.190
fourier_vec_date_period_48 lubridate_month_date_label_true_Feb
18936.14 -277734.8
s.e. 21304.62 22149.8
lubridate_month_date_label_true_Mar
-513924.41
s.e. 29139.36
lubridate_month_date_label_true_Apr
-962376.02
s.e. 36698.15
lubridate_month_date_label_true_May
-1169591.78
s.e. 43641.05
lubridate_month_date_label_true_Jun
-1197772.63
s.e. 48058.64
lubridate_month_date_label_true_Jul
-1036116.8
s.e. 48211.1
lubridate_month_date_label_true_Aug
-999032.90
s.e. 44062.88
lubridate_month_date_label_true_Sep
-1189677.89
s.e. 37261.45
lubridate_month_date_label_true_Oct
-1080367.69
s.e. 30674.67
lubridate_month_date_label_true_Nov
-817200.87
s.e. 25877.97
lubridate_month_date_label_true_Dec
-283362.29
s.e. 21328.29
sigma^2 estimated as 8.392e+09: log likelihood=-2457.31
AIC=4954.62 AICc=4959.53 BIC=5019.77Modelos para estacionalidades complejas 📈
TBATS
Una aproximación alternativa a los modelos de regresión dinámica es el modelo TBATS3334 que considera35:
- Términos trigonométricos para la estacionalidad36
- Transformaciones de Box-Cox para heterogeneidad.
- Errores ARMA para dinámica de corto plazo.
- Tendencia37
- Estacionalidad para períodos múltiples enteros y no-enteros
Show code
# TBATS
model_spec_tbats <- seasonal_reg() %>%
set_engine("tbats")
model_fit_tbats <- model_spec_tbats %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_tbats
parsnip model object
Fit time: 10.3s
TBATS(0.42, {2,1}, -, {<12,5>})
Call: forecast::tbats(y = outcome, use.parallel = use.parallel)
Parameters
Lambda: 0.419882
Alpha: 0.2311934
Gamma-1 Values: -0.001020232
Gamma-2 Values: 0.001317692
AR coefficients: 1.16936 -0.366812
MA coefficients: -0.83869
Seed States:
[,1]
[1,] 1019.049338
[2,] 97.333918
[3,] 51.268883
[4,] -3.369889
[5,] 4.590630
[6,] 4.225338
[7,] 13.800261
[8,] 16.167936
[9,] -12.308149
[10,] -1.550687
[11,] -9.391446
[12,] 0.000000
[13,] 0.000000
[14,] 0.000000
attr(,"lambda")
[1] 0.4198821
Sigma: 18.89109
AIC: 5437.111Modelo Prophet 📈
Prophet es un procedimiento de pronśtico de series de tiempo basado en modelos aditivos donde se ajustan tendencia no-lineals con estacionalidad anual, semanal y diaria, incorporando además efectos de calendario. Su potencial radica en el modelamiento de series con fuerte componente estacional y varias temporadas de datos históricos.
Prophet puede ser considerado un modelo de regresión no lineal38 de forma:
\[y_t = g(t) + s(t) + h(t) + \varepsilon_t\]
donde \(g(t)\) describe una tendencia lineal por tramos (piecewise-linear trend), \(s(t)\) describe los variados patrones estacionales, \(h(t)\) captura los efectos de calendario y \(\varepsilon_t\) es un término de ruido blanco asociado al error. Es importante destacar que:
- La cantidad de nodos (knots) o puntos de quiebre (changepoints) para la tendencia por tramos son seleccionados de manera automática si no son especificados.
- La componente estacional consiste en términos de Fourier para los períodos relevantes, usando 10 para la estacionalidad anual y 3 para la semanal.
- Los efectos de calendario son incorporados como variables dummy
- El modelo es estimado usando una aproximación Bayesiana para la selección automática de parámetros del modelo.
Show code
# PROPHET
model_spec_prophet <- prophet_reg() %>%
set_engine("prophet")
model_fit_prophet <- model_spec_prophet %>%
fit(
y ~ date,
data = training(splits)
)
model_fit_prophet
parsnip model object
Fit time: 892ms
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- changepoint.range: 0.8
- yearly.seasonality: 'auto'
- weekly.seasonality: 'auto'
- daily.seasonality: 'auto'
- seasonality.mode: 'additive'
- changepoint.prior.scale: 0.05
- seasonality.prior.scale: 10
- holidays.prior.scale: 10
- logistic_cap: NULL
- logistic_floor: NULL
- extra_regressors: 0Finalmente, adicionamos la regresión lineal pasada y el modelo estacional ingenuo como benchmarks y obtenemos:
Show code
# LM
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)
)
# SNAIVE
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: 16ms
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 2913823Show code
# MODELTIME TABLE
models_tbl <- modeltime_table(
model_fit_ses,
model_fit_holt,
model_fit_holt_damped,
model_fit_hw,
model_fit_hw_m,
model_fit_hw_m_damped,
model_fit_prophet,
model_fit_tbats,
model_fit_arima,
model_fit_linreg_arima,
model_fit_lm,
model_fit_snaive
)
models_tbl
# Modeltime Table
# A tibble: 12 × 3
.model_id .model .model_desc
<int> <list> <chr>
1 1 <fit[+]> ETS(A,N,N)
2 2 <fit[+]> ETS(A,A,N)
3 3 <fit[+]> ETS(A,AD,N)
4 4 <fit[+]> ETS(A,A,A)
5 5 <fit[+]> ETS(M,MD,M)
6 6 <fit[+]> ETS(M,MD,M)
7 7 <fit[+]> PROPHET
8 8 <fit[+]> TBATS(0.42, {2,1}, -, {<12,5>})
9 9 <fit[+]> ARIMA(1,0,2)(0,1,1)[12] WITH DRIFT
10 10 <fit[+]> REGRESSION WITH ARIMA(1,1,1)(0,0,2)[12] ERRORS
11 11 <fit[+]> LM
12 12 <fit[+]> SNAIVE [12] 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: 830 × 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 820 more rowsShow code
calibration_tbl %>%
# Grafico pronosticos
plot_modeltime_forecast(.interactive = TRUE)
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
)
Show code
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] modeltime_1.1.0 rsample_0.1.1 parsnip_0.1.7
[4] USgas_0.1.1 quantmod_0.4.18 TTR_0.24.2
[7] xts_0.12.1 zoo_1.8-9 timetk_2.6.2
[10] xaringanExtra_0.5.5 forcats_0.5.1 stringr_1.4.0
[13] dplyr_1.0.7 purrr_0.3.4 readr_2.1.1
[16] tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5
[19] tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] readxl_1.3.1 backports_1.4.0 workflows_0.2.4
[4] plyr_1.8.6 lazyeval_0.2.2 splines_4.1.2
[7] crosstalk_1.2.0 listenv_0.8.0 inline_0.3.19
[10] digest_0.6.29 foreach_1.5.1 htmltools_0.5.2
[13] yardstick_0.0.9 fansi_0.5.0 magrittr_2.0.1
[16] memoise_2.0.1 tune_0.1.6 tzdb_0.2.0
[19] recipes_0.1.17 globals_0.14.0 modelr_0.1.8
[22] gower_0.2.2 matrixStats_0.61.0 RcppParallel_5.1.4
[25] hardhat_0.1.6 forecast_8.15 tseries_0.10-49
[28] prettyunits_1.1.1 dials_0.0.10 colorspace_2.0-2
[31] rvest_1.0.2 haven_2.4.3 xfun_0.28
[34] callr_3.7.0 crayon_1.4.2 jsonlite_1.7.2
[37] survival_3.2-13 iterators_1.0.13 glue_1.5.1
[40] gtable_0.3.0 ipred_0.9-12 V8_3.6.0
[43] pkgbuild_1.2.1 rstan_2.21.2 future.apply_1.8.1
[46] scales_1.1.1 DBI_1.1.1 Rcpp_1.0.7
[49] viridisLite_0.4.0 GPfit_1.0-8 stats4_4.1.2
[52] lava_1.6.10 StanHeaders_2.21.0-7 prodlim_2019.11.13
[55] htmlwidgets_1.5.4 httr_1.4.2 ellipsis_0.3.2
[58] loo_2.4.1 pkgconfig_2.0.3 farver_2.1.0
[61] nnet_7.3-17 sass_0.4.0.9000 dbplyr_2.1.1
[64] janitor_2.1.0 utf8_1.2.2 tidyselect_1.1.1
[67] labeling_0.4.2 rlang_0.4.12 DiceDesign_1.9
[70] reactR_0.4.4 munsell_0.5.0 cellranger_1.1.0
[73] tools_4.1.2 cachem_1.0.6 cli_3.1.0
[76] generics_0.1.1 broom_0.7.10 evaluate_0.14
[79] fastmap_1.1.0 yaml_2.2.1 processx_3.5.2
[82] knitr_1.36 fs_1.5.1 future_1.23.0
[85] nlme_3.1-155 reactable_0.2.3 xml2_1.3.3
[88] compiler_4.1.2 rstudioapi_0.13 plotly_4.10.0
[91] curl_4.3.2 reprex_2.0.1 lhs_1.1.3
[94] bslib_0.3.1 stringi_1.7.6 ps_1.6.0
[97] highr_0.9 lattice_0.20-45 Matrix_1.4-0
[100] urca_1.3-0 vctrs_0.3.8 pillar_1.6.4
[103] lifecycle_1.0.1 furrr_0.2.3 lmtest_0.9-39
[106] jquerylib_0.1.4 data.table_1.14.2 R6_2.5.1
[109] gridExtra_2.3 parallelly_1.29.0 codetools_0.2-18
[112] distill_1.3 MASS_7.3-55 assertthat_0.2.1
[115] withr_2.4.3 fracdiff_1.5-1 parallel_4.1.2
[118] hms_1.1.1 quadprog_1.5-8 grid_4.1.2
[121] rpart_4.1-15 timeDate_3043.102 class_7.3-19
[124] snakecase_0.11.0 rmarkdown_2.11.3 prophet_1.0
[127] downlit_0.4.0 pROC_1.18.0 lubridate_1.8.0 feature engineering↩︎
(OLS): Ordinary Least Squares↩︎
Debemos relajar el supuesto de independencia de las muestras obtenidas a partir de la población↩︎
Este segmento es un resumen del siguiente capítulo de
fpp3: https://otexts.com/fpp3/expsmooth.html↩︎ponderadores↩︎
Una opción es definir el punto de partida como el primer valor observado, \(\ell_1 = y_1\)↩︎
De ahí el nombre del método↩︎
o importancia↩︎
No hay aprendizaje por parte del modelo↩︎
Con valores de \(\alpha\) altos, el decaimiento es lento.↩︎
Es decir, los pronósticos son una función lineal de \(h\)↩︎
Gardner & McKenzie (1985)↩︎
Asegura que los estimadores de los indices estacionales usados para pronosticar provengan del último periodo estacional de la muestra.↩︎
Es decir, hace \(m\) períodos atrás.↩︎
Además de la estimación puntual lograda con un suavizamiento exponencial.↩︎
Usualmente conocida como “All models are wrong, some are useful”↩︎
Nota del autor: ¡siempre había querido poner la cita completa en alguna parte!🤩↩︎
Aleatorio↩︎
Ya que podrían no aventajarse de la dependencia temporal en las observaciones de la manera más efectiva↩︎
También hallados en la literatura como metodología de pronóstico de Box-Jenkins↩︎
Por otro lado, recordemos que una serie de ruido blanco es estacionaria, no importa cuando se le observe, debería verse similar en cualquier punto del tiempo↩︎
La ACF de una serie de tiempo estacionaria decae a cero relativamente rápido, mientras que para series no estacionarias lo hace lentamente↩︎
¡Pronóstico ingenuo!↩︎
Hablando de la diferencia usual como primera diferencia↩︎
Un test útil sería el test de Kwiatkowski-Phillips-Schmidt-Shin (KPSS)↩︎
Si la media de no es cero, entonces se reemplaza \(y_t\) por \(y_t - \mu\), donde \(\mu\) es la media de la serie, en este caso \(c\)↩︎
No estacional↩︎
incluso más de una vez↩︎
Paper: Automatic time series forecasting: the forecast package for R↩︎
En este caso, en un año.↩︎
Feriados, información de la competencia, etc.↩︎
T: estacionalidad trigonométrica, B: transformación Box-Cox, A: errores ARIMA, T: tendencia, S: componentes estacionales↩︎
Paper: Forecasting Time Series With Complex Seasonal Patterns Using Exponential Smoothing↩︎
Puedes encontrar una buena comparación del modelo TBATS con regresión Dinámica acá: Forecasting weekly data↩︎
Admite múltiples estacionalidades↩︎
Usualmente amortiguada↩︎
Similar a los Modelos Aditivos Generalizados (GAM)↩︎