07 Time Series Models
Eduardo Villarreal
13/8/2021
require(tidyverse)
require(data.table)
require(gridExtra) #Libreria para hacer paneles de gráficas con ggplot2
require(forecast) #Libreria para análisis de series de tiempo
require(fpp2) #Datos de eries de tiempo
require(timetk) #To work with time series
require(modeltime) #To model time series models
require(tidyquant) #Timeserioes Operarions
require(hts) #Library for multiple time series data
require(zoo)En esta ocasión vamos a utilizar datos de CoVid19 para predecir la tendencia de nuevos casos de la serie diaria. La fuente de información es Our World in Data
#Descargar los datos de CoVid
covid_wdb = fread('https://covid.ourworldindata.org/data/owid-covid-data.csv')
#crear el campo de Fecha
covid_wdb$Date = as.Date(covid_wdb$date, format = '%Y-%m-%d')
head(covid_wdb)1 Análisis descriptivo de la serie
La variable a la cual haremos referencia es new_cases_per_million y estamos interesados en el total de casos para México, Colombia, Peru y Brasil
pais = c('Mexico', 'Colombia', 'Peru', 'Brazil')
pais_MA = paste(pais, 'MA14', sep = '_')
#Querie para filtrar solo los paises de interés
MA_covid = covid_wdb %>%
filter(location %in% pais) %>%
mutate(new_cases = new_cases_per_million) %>%
dplyr::select(Date, location, new_cases) %>%
spread(location, new_cases, fill = 0)
#Crear Total y solo tomar desde Abril del 2020
MA_covid$Total = rowSums(MA_covid[, 2 : ncol(MA_covid)])
MA_covid = MA_covid %>%
filter(Date >= '2020-04-01') %>%
gather(location, new_cases, 2 : ncol(.))
#Grafica de los casos diarios
MA_covid %>%
ggplot(aes(x = Date, y = new_cases, color = location)) +
geom_line() +
facet_wrap(~location, scales = 'free', nrow = 1) +
theme_bw() +
scale_color_viridis_d(begin = 0.2, end = 0.8) +
labs(title = 'Covid New cases/million by location') +
theme(legend.position = 'bottom')
1.1 Seasonal Decomposition
Como vimos en la introducción, una serie de tiempos es la combinación de 3 elementos:
Nivel
Tendencia
Estacionalidad
\(y_t = N_t + T_t + S_t +\epsilon_t\)
Ajuste por Estacionalidad
Si removemos el componente estacional de los datos el resultado es una serie desestacionalida:
\(y_t - S_t = N_t + T_t + \epsilon_t\) si la series es aditiva y si la series es múltiplicativa, obtenemos:
\(\frac{y_t}{S_t} = N_t \times T_t \times \epsilon_t\)
Desestacionalizar una serie es útil cuando la variación estacional de la serie no es de interés. Por ejemplo, la tasa mensual de desempleo suele analizarse sin estacionalidad. Por ejemplo, una caÃda en la tasa de desempleo debido a los trabajos escolares.
Siguiendo con la serie de Covid Total, para ajustar por estacionalidad primero aplicamos la función decompose y posteriormente aplicamos la función seasadj:
MA_covid = MA_covid %>%
spread(location, new_cases, fill = 0)
covid_total = ts(MA_covid$Total, frequency = 7, start = c(2020, 4))
covid_dcomposed = decompose(covid_total)
autoplot(covid_dcomposed) +
theme_bw()
El objeto covid_decomsed tiene los componentes de Tendencia, Etacionalidad y el componente Aleatorio. Por ejemplo, para mejorar la visualización de la gráfica, podemos eliminar el ciclo estacional haciendo:
covid_dcomposed = seasadj(covid_dcomposed)
autoplot(covid_dcomposed) +
theme_bw() +
labs(title = 'Serie Covid Total Desestacionalizada')
En este caso estamos observando la serie \(z_t = y_t - S_t\); es decir, el componente de tendencia y error aleatorio.
Si quisiéramos observar la tendencia de la serie, podemos extraer la tendencia como:
covid_dcomposed = decompose(covid_total)
covid_trend = covid_dcomposed$trend
autoplot(covid_trend) +
theme_bw() +
labs(title = 'Tendencia Serie Covid Total')
1.2 Moving Average Smoothing
Otra forma de visualizar la serie es utilizando Moving Averages.
Este tipo de computos se refiere a ventanas de tiempo bajo las cuales se calcula algún tipo de estadÃstico como la media la desviación estándar o la correlación entre dos series de tiempo. Un ejemplo es el cómputo de Media Móvil o MA(L):
\(MA(L) = \frac{1}{L}\sum_{t=1}^L X_t\)
MA_smooth = MA_covid %>%
select(Date, Total) %>%
tq_mutate(
select = Total,
mutate_fun = rollapply,
width = 20,
FUN = mean,
col_rename = 'MA20') %>%
tq_mutate(
select = Total,
mutate_fun = rollapply,
width = 50,
FUN = mean,
col_rename = 'MA50')## Warning: `type_convert()` only converts columns of type 'character'.
## - `df` has no columns of type 'character'
## Warning: `type_convert()` only converts columns of type 'character'.
## - `df` has no columns of type 'character'#Graficar
MA_smooth %>%
ggplot(aes(x = Date, y = Total)) +
geom_line(alpha = 0.6, col = 'grey50') +
geom_line(aes(y = MA20), size = 1, col = 'blue') +
geom_line(aes(y = MA50), size = 1, col = 'red') +
theme_bw() +
labs(title = 'Suavisamiento con MA20 y MA50')## Warning: Removed 19 row(s) containing missing values (geom_path).## Warning: Removed 49 row(s) containing missing values (geom_path).
En este caso, estamos usando un MA(20) como una medida de la tendencia a corto plazo y un MA(50) como una medida de largo plazo. Lo interesante de estos patrones son los cruces entre las dos lÃneas:
Cuando MA(20) > MA(50); la tendencia de casos es a la alza
Cuando MA(20) > MA(50); la tendencia de casos es a la baja
cuando MA(20) = MA(50) hay un cambio de tendencia
Lo que podemos ver acá es que nuestro MA(50) empieza a tener un momentum negativo acelerado, indicando una posible 5ta ola de contagios en la región de interés.
1.3 Caso particular de México
Ahora vamos a hacer el caso particular de México
MA_smooth = MA_covid %>%
select(Date, Mexico) %>%
tq_mutate(
select = Mexico,
mutate_fun = rollapply,
width = 21,
FUN = mean,
col_rename = 'MA_short') %>%
tq_mutate(
select = Mexico,
mutate_fun = rollapply,
width = 21 * 3,
FUN = mean,
col_rename = 'MA_long')## Warning: `type_convert()` only converts columns of type 'character'.
## - `df` has no columns of type 'character'
## Warning: `type_convert()` only converts columns of type 'character'.
## - `df` has no columns of type 'character'#Graficar
MA_smooth %>%
ggplot(aes(x = Date, y = Mexico)) +
geom_line(alpha = 0.6, col = 'grey50') +
geom_line(aes(y = MA_short), size = 1, col = 'blue') +
geom_line(aes(y = MA_long), size = 1, col = 'red') +
theme_bw() +
labs(title = 'Suavisamiento con MA_short y MA_long')## Warning: Removed 20 row(s) containing missing values (geom_path).## Warning: Removed 62 row(s) containing missing values (geom_path).
Qué podemos observar del comportamiento en la 3ra ola?
1.4 STL Decomposition
STL es un acrónimo de Seasonal and Trend decomposition using LOESS en dondo LOESS es un método para estimar relaciones no lineales que usa regression splins para ajustarse en forma adaptativa.
STL tiene muchas ventajas:
Puede manejar cualquier tipo de estacionalidad: diaria, semanal, anual, etc.
El componente estacional es adaptativo en el tiempo
Es robusto ante la presencia de outliers
Hay dos parámetros que debe poner el usiario: t.window (trend-cycle window) y s.window (seasonal window). Estos dos parámetros controlan qué tan rápido pueden adaptarse ambos componentes. Valores pequeños permiten cambios más rápidos. Ambos valores deben ser impares.
Una interpretación más precisa de ambos parámetros es t.windowes el número de observaciones consecutivas utilizadas para estimar el componente de ciclo-tendencia.
MA_covid$Total %>%
ts(frequency = 7) %>%
mstl() %>%
autoplot()
Observe cómo, a diferencia de decompose, el componente estacional varia a lo largo del tiempo.
Si quisieramos ajustar los datos eliminando la estacionalidad, podemos hacer \(z_t = y_t - S_t\):
seasonal_eff = MA_covid$Total %>%
ts(frequency = 7) %>%
mstl()
#Seasonal djusted series
seas_adj = MA_covid$Total - seasonal_eff[, 3]
autoplot(seas_adj) +
theme_bw() +
labs(title = 'Seasonal adj Series using STL')
Y si quisieramos observar la tendencia, podemos hacer:
seas_adj = MA_covid$Total - seasonal_eff[, 3] - seasonal_eff[, 4]
autoplot(seas_adj) +
theme_bw() +
labs(title = 'Trend Component Series using STL')
1.4.1 Qué tan fuerte es la Tendencia y la Estacionalidad?
Para datos con una tendencia muy fuerte los datos ajustados por estacionalidad deben tener mucha más varianza que el componente aleatorio. AsÃ, el ratio \(\text{Var}(\epsilon_t) / \text{Var}(T_t + \epsilon_t)\) deberÃa ser relativamente pequeño. Para datos con poca o nula tendencia ambas variazas deberian ser aproximadamente iguales. Entonces definimos la fuerza de la tendencia \(F_t\) como:
\(F_T = max(0, 1 - \frac{\text{Var}(\epsilon_t)}{\text{Var}(T_t + \epsilon_t)})\)
De manera similar podemos definir la fortaleza del factor estacional como:
\(F_S = max(0, 1 - \frac{\text{Var}(\epsilon_t)}{\text{Var}(S_t + \epsilon_t)})\)
#Funcion para definir la fortaleza de la Tendencia y la Estacionalidad
FT_strenght = function(data, freq = 12){
x = data
x = ts(x, frequency = freq)
#Time Series decomposition using STL
x_decompose = mstl(x)
trend_t = x_decompose[, 2]
season_t = x_decompose[, 3]
reminder_t = x_decompose[, 4]
#Strength of a Trend
F_Trend = var(reminder_t) / var(trend_t + reminder_t)
F_Trend = max(0, 1 - F_Trend)
return(F_Trend)
}
FS_strenght = function(data, freq = 12){
x = data
x = ts(x, frequency = freq)
#Time Series decomposition using STL
x_decompose = mstl(x)
trend_t = x_decompose[, 2]
season_t = x_decompose[, 3]
reminder_t = x_decompose[, 4]
#Strength of a Seasonal
F_Season = var(reminder_t) / var(season_t + reminder_t)
F_Season = max(0, 1 - F_Season)
return(F_Season)
}
#Caso particular de Covid total
FS_strenght(MA_covid$Total, freq = 7)## [1] 0.5368644FT_strenght(MA_covid$Total, freq = 7)## [1] 0.8571836#Caso de todas las series
sapply(MA_covid[, 2 : ncol(MA_covid)], FS_strenght)## Brazil Colombia Mexico Peru Total
## 0.06579056 0.11419301 0.09722337 0.18998776 0.12207984sapply(MA_covid[, 2 : ncol(MA_covid)], FT_strenght)## Brazil Colombia Mexico Peru Total
## 0.6137798 0.9669521 0.6750007 0.3625613 0.7518501Esto es muy útil cuando queremos detectar en un conjunto de series de tiempo quiénes tienen Tendecia y Estacionalidad fuerte.
1.5 Descomposición Estacional con Series de Fourier
Una serie de Fourier es una serie infinita que converge puntualmente a una función periódica y continua a trozos (o por partes). Las series de Fourier constituyen la herramienta matemática básica del análisis de Fourier empleado para analizar funciones periódicas a través de la descomposición de dicha función en una suma infinita de funciones sinusoidales mucho más simples (como combinación de senos y cosenos con frecuencias enteras).
El nombre se debe al matemático francés Jean-Baptiste Joseph Fourier, que desarrolló la teorÃa cuando estudiaba la ecuación del calor. Fue el primero que estudió tales series sistemáticamente, y publicó sus resultados iniciales en 1807 y 1811. Esta área de investigación se llama algunas veces análisis armónico.
Las series de Fourier tienen la forma:
\(y_t = \sum_{k=1}^K \left[\gamma_k \sin\left(\textstyle\frac{2\pi t k}{m}\right) \psi_k\cos\left(\textstyle\frac{2\pi t k}{m}\right)\right]\)
En donde \(1 < K < m/2\). Mientras más grande es el valor de \(K\) el componente estacional es más complejo. El valor óptimo de \(K\) puede estimase minimizando algún criterio como AIC o leave-one-out CV
La función fourier() puede usarse para estimar la serie de fourier:
#Estimación de la serie de fourier
Fourier = MA_covid$Total %>%
ts(frequency = 7) %>%
fourier(K = 2) %>%
as.data.frame()
#Estimacion de los coeficientes de Fourier usando regresión lineal
temp = Fourier
temp$y = MA_covid$Total
seasonal_model = lm(y ~., data = temp)
#Seasonal component
temp$seasonal = seasonal_model$fitted.values
#graficar el componente esacional
plot(temp$seasonal, type = 'l')
1.5.1 Descomposición con Regresión Lineal y Series de Fourier
Podemos ahora descomponer la serie como:
\(y_t = \alpha + \beta t + \sum_{k=1}^K \left[\gamma_k \sin\left(\textstyle\frac{2\pi t k}{m}\right) \psi_k\cos\left(\textstyle\frac{2\pi t k}{m}\right)\right] + \varepsilon_t\)
#Estimat con Regresión lineal
temp = MA_covid %>%
select(Date, Total) %>%
mutate(Trend = 1 : nrow(.))
#Estimacion de Series de Fourier
Fourier = MA_covid$Total %>%
ts(frequency = 7) %>%
fourier(K = 2) %>%
as.data.frame()
#Aumentar el data frame
temp = bind_cols(temp, Fourier)
head(temp)#Modelo de Regresión
model_decompose = lm(Total ~., data = temp[, -1])
summary(model_decompose)##
## Call:
## lm(formula = Total ~ ., data = temp[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -623.68 -128.13 -13.89 126.19 732.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 237.58659 17.99823 13.201 < 2e-16 ***
## Trend 1.27250 0.06201 20.522 < 2e-16 ***
## `S1-7` 95.94832 12.69790 7.556 2.02e-13 ***
## `C1-7` -45.00394 12.71754 -3.539 0.00044 ***
## `S2-7` 29.90210 12.70626 2.353 0.01900 *
## `C2-7` 28.05268 12.70873 2.207 0.02775 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 201.3 on 496 degrees of freedom
## Multiple R-squared: 0.5018, Adjusted R-squared: 0.4968
## F-statistic: 99.93 on 5 and 496 DF, p-value: < 2.2e-16Ahora podemos descomponer la serie:
trend = coef(model_decompose)[1] + coef(model_decompose)['Trend'] * temp$Trend
season = temp$Total - trend - residuals(model_decompose)
components = cbind(
time = temp$Trend,
data = temp$Total,
trend = trend,
season = season,
reminder = residuals(model_decompose))
#graficar
components %>%
as.data.frame() %>%
gather(Series, Value, 2 : ncol(.)) %>%
ggplot(aes(x = time, y = Value)) +
geom_line() +
facet_wrap(~Series, scales = 'free', ncol = 1)
Podemos tambien mejorar la tendencia utilizando splines
Podemos particionar el tiempo (la tendencia) en 3 momentos diferentes, definidos arbitrariamente como \(Q_1, Q_2, Q_3\) que representan el quartil 1, 2 y 3 respectivamente:
summary(temp$Trend)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 126.2 251.5 251.5 376.8 502.0Nuestra partición quearÃa como c(126, 251, 376) y usamos el concepto de Cubic Splin para ajustar la tendencia:
temp = MA_covid %>%
select(Date, Total) %>%
mutate(Trend = 1 : nrow(.))
#compute Fourier Series
Fourier = MA_covid$Total %>%
ts(frequency = 7) %>%
fourier(K = 2) %>%
as.data.frame()
temp = bind_cols(temp, Fourier)
require(splines)## Loading required package: splines#compute Trend using basis functions
m2 = lm(Total ~ bs(Trend, knots = c(126, 251,376)) - 1, data = temp)
Trend = m2$fitted.values
summary(m2)##
## Call:
## lm(formula = Total ~ bs(Trend, knots = c(126, 251, 376)) - 1,
## data = temp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -540.98 -108.18 -7.81 88.13 830.47
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## bs(Trend, knots = c(126, 251, 376))1 204.45 44.93 4.551 6.73e-06 ***
## bs(Trend, knots = c(126, 251, 376))2 709.20 49.40 14.356 < 2e-16 ***
## bs(Trend, knots = c(126, 251, 376))3 329.63 43.43 7.591 1.59e-13 ***
## bs(Trend, knots = c(126, 251, 376))4 891.88 51.98 17.157 < 2e-16 ***
## bs(Trend, knots = c(126, 251, 376))5 1050.69 54.18 19.393 < 2e-16 ***
## bs(Trend, knots = c(126, 251, 376))6 318.71 49.37 6.456 2.57e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 174.5 on 496 degrees of freedom
## Multiple R-squared: 0.9232, Adjusted R-squared: 0.9222
## F-statistic: 993.2 on 6 and 496 DF, p-value: < 2.2e-16#Decompose the time series:
temp$Trend = Trend
temp$seq = 1 : nrow(temp)
model_decompose2 = lm(Total ~., data = temp[, -1])
trend = coef(model_decompose2)[1] + coef(model_decompose2)['Trend'] * temp$Trend
season = temp$Total - trend - residuals(model_decompose2)
components = cbind(
time = temp$seq,
data = temp$Total,
trend = trend,
season = season,
reminder = residuals(model_decompose2))
#graficar
components %>%
as.data.frame() %>%
gather(Series, Value, 2 : ncol(.)) %>%
ggplot(aes(x = time, y = Value)) +
geom_line() +
facet_wrap(~Series, scales = 'free', ncol = 1)
2 Forecast con Regresión Lineal
El enfoque más directo para modelar y proectar series de tiempo es la regresión lineal. Podemos escribir una series de tiempos \(y_t\) como:
\(y_t = \beta_0 + \beta_1 T_t + \sum_{s=1}^S \phi_s S_t +\varepsilon_t\)
En donde \(S_t\) es el componente estacional que en el caso más sencillo es una variabel dummy. Para más detalles sobre regresión y variabels dummy pudes revisar https://rpubs.com/blad0914.
El paquete forecast trae consigo la función tslm que acepta una formula como objeto y estima los coeficientes de regresión.
covid_total = MA_covid$Total %>%
ts(frequency = 7)
m_lm = tslm(covid_total ~ trend + season)
summary(m_lm)##
## Call:
## tslm(formula = covid_total ~ trend + season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -644.8 -131.7 -10.8 122.1 711.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 297.99640 28.29678 10.531 < 2e-16 ***
## trend 1.27207 0.06195 20.535 < 2e-16 ***
## season2 5.14169 33.52103 0.153 0.8782
## season3 25.31181 33.52120 0.755 0.4505
## season4 -37.44749 33.52149 -1.117 0.2645
## season5 -134.93887 33.52189 -4.025 6.58e-05 ***
## season6 -220.31914 33.63893 -6.550 1.45e-10 ***
## season7 -59.91021 33.63915 -1.781 0.0755 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 201.1 on 494 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.4978
## F-statistic: 71.94 on 7 and 494 DF, p-value: < 2.2e-16Es importate notar que trend y season son creados en automático por la función tslm. Con el modelo ya calibrado, ahora, podemos hacer un forecast con la función furecast para los siguientes 12 dias:
fcs_lm = forecast(m_lm, h = 12)
fcs_lm## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72.71429 717.5280 456.8476 978.2084 318.4071 1116.649
## 72.85714 879.2090 618.5286 1139.8894 480.0881 1278.330
## 73.00000 940.3913 679.7144 1201.0682 541.2756 1339.507
## 73.14286 946.8051 686.1281 1207.4820 547.6894 1345.921
## 73.28571 968.2472 707.5703 1228.9241 569.1316 1367.363
## 73.42857 906.7600 646.0831 1167.4369 507.6444 1305.876
## 73.57143 810.5407 549.8638 1071.2176 411.4250 1209.656
## 73.71429 726.4325 465.7088 987.1562 327.2452 1125.620
## 73.85714 888.1135 627.3898 1148.8372 488.9262 1287.301
## 74.00000 949.2958 688.5749 1210.0166 550.1128 1348.479
## 74.14286 955.7095 694.9887 1216.4304 556.5266 1354.892
## 74.28571 977.1517 716.4309 1237.8726 577.9688 1376.335autoplot(fcs_lm) +
theme_bw()
Y también podemos computar los métricos de frecast accuracy con la función accuracy:
accuracy(m_lm)## ME RMSE MAE MPE MAPE MASE
## Training set -1.439895e-14 199.5168 155.4599 -45.38192 65.56893 1.236899
## ACF1
## Training set 0.5293117Para llevar un registro de los distintos pronósticos que vamos a hacer, vamos a crear un Data Frame de resultados con la fecha, los datos actuales y los estimados:
FCS_Results = MA_covid %>%
select(Date, Total) %>%
mutate(m_lm = m_lm$fitted.values)Un aspecto importante del forecast es que hay algunos supuestos que deben cumplir los residuales:
los residuales \(\varepsilon_t\) deben ser \(\text{NID}(0, \sigma^2)\)
los residuales no deben estar autocorrelacionados
Para revisar estos dos spuestos podemos hacer:
checkresiduals(m_lm)
##
## Breusch-Godfrey test for serial correlation of order up to 14
##
## data: Residuals from Linear regression model
## LM test = 292.09, df = 14, p-value < 2.2e-16Aún cuando la distribución de los residuales parece seguir una forma gausiana vemos que en el tiempo no siguen un patrón de ruido blanco y la función de autocorrelación o ACF muestra que existe autocorrelación significativa.
Este es un problema común cuando queremos modelar series de tiempo con regresión lineal.
3 Forecast con Descomposición STL
Una de las ventajas de la descomposición STL es que podemos hacer un forecast de una manera relativamente directa:
m_stl = MA_covid$Total %>%
ts(frequency = 7) %>%
mstl()
#Hacer el forecast
fcs_stl = forecast(m_stl, h = 12)
autoplot(fcs_stl)
#Revisar los residuales
checkresiduals(fcs_stl)## Warning in checkresiduals(fcs_stl): The fitted degrees of freedom is based on
## the model used for the seasonally adjusted data.
##
## Ljung-Box test
##
## data: Residuals from STL + ETS(M,A,N)
## Q* = 173.54, df = 10, p-value < 2.2e-16
##
## Model df: 4. Total lags used: 14#Gardar el valor ajustado
FCS_Results$m_stl = fcs_stl$fitted
#Computar accuracy
accuracy(fcs_stl)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -3.381969 117.7983 82.00973 -4.82153 16.98719 0.6525008 -0.20815654 Forecasting con Suavisamiento Exponencial
El Suavisamiento Exponencial es una técnica utilizada desde 1957 creada por Holt, Brown y Winters.
4.1 Suavisamiento Exponencial Simple (SES)
El caso más sencillo es el suavisamiento exponencial simple. Este método es adecuado cuando no hay tendencia y estacionalidad.
El forecast en el tiempo \(T+1\) puede escribirse como:
\(\hat{y}_{T+1|t} = \alpha y_T + (1-\alpha) \hat{y}_{T|T-1}\)
en donde el parámetro de suavisamiento \(0 < \alpha < 1\)
Si \(\ell_{t}\) es el nivel de la serie en el tiempo \(t\), entonces podemos escribir el modelo en forma estructural como:
\(\text{Forecast: } \hat{y}_{t+h|t} = \ell_{t}\)
\(\text{Smoothing equation: } \ell_{t} = \alpha y_{t} + (1 - \alpha)\ell_{t-1}\)
4.1.1 Estimación del modelo
Para estimar el modelo son necesarion métodos de óptimización no lineal con restricciones. Pero en R podemos hacerlo con la libreria forecast a través de la función ses:
m_ses = MA_covid$Total %>%
ts(frequency = 7) %>%
ses(h = 12)
summary(m_ses)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = ., h = 12)
##
## Smoothing parameters:
## alpha = 0.1598
##
## Initial states:
## l = 27.5591
##
## sigma: 165.7849
##
## AIC AICc BIC
## 8256.867 8256.915 8269.523
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4.715914 165.4543 120.9838 -5.948317 25.62998 0.9625934 0.05639507
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72.71429 405.9303 193.4684 618.3921 80.99792 730.8626
## 72.85714 405.9303 190.7719 621.0886 76.87398 734.9865
## 73.00000 405.9303 188.1088 623.7517 72.80108 739.0594
## 73.14286 405.9303 185.4778 626.3827 68.77739 743.0831
## 73.28571 405.9303 182.8779 628.9826 64.80115 747.0594
## 73.42857 405.9303 180.3080 631.5526 60.87073 750.9898
## 73.57143 405.9303 177.7669 634.0936 56.98457 754.8759
## 73.71429 405.9303 175.2539 636.6066 53.14123 758.7193
## 73.85714 405.9303 172.7680 639.0925 49.33930 762.5212
## 74.00000 405.9303 170.3082 641.5523 45.57748 766.2830
## 74.14286 405.9303 167.8739 643.9866 41.85453 770.0060
## 74.28571 405.9303 165.4643 646.3962 38.16927 773.6912En este caso, el problema de optimización tiene que lidear con 2 valores para estimar:
El valor del parámetro de suavizamiento, que en este caso es \(\alpha = 0.1576\)
El valor inicial \(\ell_0 = 19.90\)
4.1.2 Forecasting con el modelo SES
Una vez estimado el modelo, podemos hacer el pronóstico:
fcs_ses = forecast(m_ses)
autoplot(fcs_ses)
Es importante que por definición matemática el valor pronósticado es y será siempre constante.
Ahora podemos registrar el valor estimado y analizar los reidales
FCS_Results$m_ses = fcs_ses$fitted
checkresiduals(m_ses)
##
## Ljung-Box test
##
## data: Residuals from Simple exponential smoothing
## Q* = 367.42, df = 12, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 144.2 Suavisamiento de Holt
Podemos ahora incluir el componente de tendencia:
\(\text{Forecast equation: } \hat{y}_{t+h|t} = \ell_{t} + hb_{t}\)
\(\text{Level equation: }\ell_{t} = \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\)
\(\text{Trend equation: }b_{t} = \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1}\)
Para asegurar estacionariedad los valores de \(0 < \alpha < 1\) y \(0 < \beta < 1\)
Igual que el caso anterior, podemos estimar y hacer el forecast para los siguientes 12 dias:
m_holt = MA_covid$Total %>%
ts(frequency = 7) %>%
holt(h = 12)
summary(m_holt)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = ., h = 12)
##
## Smoothing parameters:
## alpha = 0.1249
## beta = 0.0155
##
## Initial states:
## l = 6.0809
## b = 2.4088
##
## sigma: 166.0007
##
## AIC AICc BIC
## 8260.161 8260.282 8281.254
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1.424795 165.3381 122.5772 -7.353286 25.22181 0.9752711
## ACF1
## Training set 0.07380583
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72.71429 349.0439 136.305348 561.7824 23.6883808 674.3993
## 72.85714 340.3918 125.567114 555.2166 11.8457700 668.9379
## 73.00000 331.7398 114.371504 549.1082 -0.6963378 664.1760
## 73.14286 323.0878 102.685244 543.4904 -13.9888297 660.1645
## 73.28571 314.4358 90.479948 538.3917 -28.0751188 656.9467
## 73.42857 305.7838 77.732417 533.8352 -42.9906846 654.5583
## 73.57143 297.1318 64.424770 529.8388 -58.7628742 653.0264
## 73.71429 288.4798 50.544398 526.4151 -75.4109710 652.3705
## 73.85714 279.8277 36.083752 523.5717 -92.9465197 652.6020
## 74.00000 271.1757 21.039990 521.3115 -111.3738694 653.7253
## 74.14286 262.5237 5.414503 519.6329 -130.6908887 655.7383
## 74.28571 253.8717 -10.787619 518.5310 -150.8897953 658.6332Ahora notemos que existen 4 estimaciones: \(\alpha = 0.113\), \(\beta = 0.021\) y \(\ell_0 = 5.93, b_0 = 2.07\).
El forecast y el análisis de los residuales es:
fcs_holt = forecast(m_holt)
FCS_Results$m_holt = fcs_holt$fitted
autoplot(fcs_holt)
checkresiduals(m_holt)
##
## Ljung-Box test
##
## data: Residuals from Holt's method
## Q* = 362.62, df = 10, p-value < 2.2e-16
##
## Model df: 4. Total lags used: 145 Suavisamiento Exponencial de Holt-Winters
Ahora vamos a incluir el tercer componente \(S_t\):
\(\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}\)
De la misma forma, los parámetros a estimar del modelo están restringidos para asegurar estacionariedad. Para hacer la estimacón del modelo usamos la función hw:
m_hw = MA_covid$Total %>%
ts(frequency = 7) %>%
hw(h = 12, seasonal = 'additive')
summary(m_hw)##
## Forecast method: Holt-Winters' additive method
##
## Model Information:
## Holt-Winters' additive method
##
## Call:
## hw(y = ., h = 12, seasonal = "additive")
##
## Smoothing parameters:
## alpha = 0.1356
## beta = 0.0183
## gamma = 0.0804
##
## Initial states:
## l = -40.6469
## b = 3.8737
## s = -27.7059 -125.8333 -47.9088 -7.3137 58.5825 76.6019
## 73.5774
##
## sigma: 137.6683
##
## AIC AICc BIC
## 8079.162 8079.800 8129.785
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.523458 136.1516 96.96892 -3.0665 26.22563 0.7715219 -0.1655926
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72.71429 104.47936 -71.94969 280.9084 -165.345580 374.3043
## 72.85714 328.47387 149.96898 506.9788 55.474197 601.4735
## 73.00000 413.83303 232.76356 594.9025 136.911180 690.7549
## 73.14286 407.69518 223.53642 591.8539 126.048669 689.3417
## 73.28571 409.14011 221.33792 596.9423 121.921451 696.3588
## 73.42857 348.68642 156.66410 540.7087 55.013620 642.3592
## 73.57143 157.28098 -39.55387 354.1158 -143.751939 458.3139
## 73.71429 33.79476 -172.17991 239.7694 -281.216316 348.8058
## 73.85714 257.78927 45.90197 469.6766 -66.264386 581.8429
## 74.00000 343.14843 124.75041 561.5464 9.137484 677.1594
## 74.14286 337.01058 111.50948 562.5117 -7.863582 681.8847
## 74.28571 338.45551 105.26857 571.6424 -18.173142 695.0842El forecast para los siguientes 12 dias es:
fcs_hw = forecast(m_hw)
autoplot(fcs_hw)
FCS_Results$m_hw = fcs_hw$fitted
checkresiduals(fcs_hw)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' additive method
## Q* = 84.744, df = 3, p-value < 2.2e-16
##
## Model df: 11. Total lags used: 146 Recapitulación de los métodos revisados
Ahora es un bue momento para recapitular y ver cómo los modelos se comportan. Para esto vamos a hacer una gráfica del Actual y los valores ajustados:
FCS_Results %>%
ggplot(aes(x = Date, y = Total)) +
geom_line(color = 'grey80') +
geom_line(aes(y = m_lm), color = 'red', lty = 2) +
geom_line(aes(y = m_stl), color = 'darkred') +
geom_line(aes(y = m_ses), color = 'darkorange') +
geom_line(aes(y = m_holt), color = 'darkblue') +
geom_line(aes(y = m_hw), color = 'green') +
theme_bw()
Y ahora vamos a imprimir el accuracy de cada uno de ellos:
fcs_acc = rbind(accuracy(m_lm),
accuracy(m_ses),
accuracy(m_holt),
accuracy(m_hw))
#Convert to Dataframe
fcs_acc = as.data.frame(fcs_acc)
fcs_acc$model = c('lm', 'ses', 'holt', 'hw')
fcs_accCual modelo utilizarÃas para hacer el pronóstico?
7 El engaño de tomar el 100% de los datos
En realidad, para que la competencia sea justa y equilibrada, los modelos deberÃan competir con datos cigos que nunca han visto.
Ahora vamos a poner a competir el modelo de holt-winter con un MAPE = 32% contra el modelo de holt que tiene un MAPE = 28% que en teorÃa es el que menor error absoluto porcentual tiene. Para esto, vamos a dejar el 95% de los datos para entrenar el modelo y el 5%% para validarlo:
#Crear el training (476 obs) y test set (25 obs)
full_set = ts(MA_covid$Total, frequency = 7)
train_set = subset(full_set, end = length(full_set) - 25)
test_set = subset(full_set, start = length(full_set) - 24)
#Computar los modelos con el training set
#Holt
m_holt = train_set %>%
holt(h = 25)
m_holt %>%
forecast(h = 25) %>%
autoplot() + autolayer(test_set) +
theme_bw()
#Holt-Winters
m_hw = train_set %>%
hw(h = 25)
m_hw %>%
forecast(h = 25) %>%
autoplot() + autolayer(test_set) +
theme_bw()
#Accuracy en Test Set para holt
accuracy(m_holt, test = test_set)## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -71.46494 213.0437 174.7907 -18.36175 30.23443 0.4171652 0.2178601#Accuracy en Test Set para HW
accuracy(m_hw, test = test_set)## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -32.29922 157.7501 126.7244 -9.896882 20.46035 0.3024476 0.137925Ahora lo que observamos es que el Holt-Winter con un 34% en el test set supera al modelo de Holt el cual tiene un 47% de error.