Pronostico con regresión lineal
Daniela Troncoso, Karla Orejarena, Miguel Cabrera
2022-09-09
El presente trabajo es un resumen detallado de los aspectos tratados en el texto Hands On Time Series Analysis with R en el que se replicaron los códigos para el aprendizaje de la estadística implicada. Este informe fue realizado por: Daniela Troncoso, Karla Orejarena y Miguel Cabrera.
¿Qué es la regresión lineal?
La regresión lineal es un modelo que pretende cuantificar la relación lineal que hay entre una variable dependiente, \(y\), y una variable o variables independientes, \(x\). Básicamente, el modelo expresa la variable dependiente como una combinación lineal de las variables independientes. Entonces, una relación lineal entre las variables dependientes e independientes puede generalizarse por medio de las siguientes ecuaciones:
En el caso de una sola variable independiente, la ecuación es la siguiente: \[Y_i= \beta_0+\beta_1 \cdot X_{1,i}+ \in_i \] Para n variables independientes, la ecuación se ve de la siguiente manera: \[Y_i =\beta_0 +\beta_1 X_{1,i}+\beta_2X_{2,i}+...+\beta_nX_{n,i}+\in_i \] Donde:
- \(i\) representa el índice de la observación \(i = 1,..., N\)
- \(Y_i\) representa la estimación de la variable dependiente \(i\)
- \(X_{j,i}\) representa el valor i de la variable independiente \(j\), donde \(j = 1,..., n\)
- \(\beta_0\) representa la estimación del término constante (o intercepción)
- \(\beta_j\) son la estimación de los parámetros (o coeficientes) correspondientes de las n variables independientes
- \(\in\) define el término de error, que no es más que toda la información que no se capturado por variables independientes para la \(i\) observación
La estimación de los coeficientes del modelo se basa en los siguientes dos pasos:
- Definir una función de costo (también conocida como función de pérdida)
- Aplicar optimización matemática para minimizar la función de costo
El enfoque de estimación más común es aplicar el método de mínimos cuadrados ordinarios (OLS) como una técnica de optimización para minimizar la suma residual de cuadrados (\(\sum_{i=1}^{N}\in_i^2\)).
El cuadrado de los residuos tiene dos efectos:
- Evita que los valores positivos y negativos se cancelen entre sí al sumarlos
- El efecto del cuadrado proporciona una penalización exponencial para los residuos con mayor distancia, ya que su costo se vuelve más alto.
Estimación de coeficientes con el método OLS
Características de la ingeniería de los componentes de la serie
Antes de crear las entradas de regresión que representan la tendencia de la serie y los componentes estacionales, primero tenemos que entender su estructura. Por lo tanto, en el libro se nos presentan ejemplos en los que se demuestra el proceso de creación de nuevas características a partir de series existentes tomadas de una base de datos llamada Usgas.
Una vez que la serie se transformó en una estructura aditiva, la transformación a una formación de regresión lineal es sencilla y sigue el mismo proceso descrito anteriormente:
library(TSstudio)
data("USgas")
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Además, se han representado las principales características de la serie utilizando la función ts_info:
ts_info(USgas)## The USgas series is a ts object with 1 variable and 238 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2019 10Como se puede ver en la trama de la serie USgas es una serie mensual con un fuerte componente estacional mensual y una línea de tendencia bastante estable. Podemos explorar la estructura de los componentes de la serie con la función ts_decompose más adelante:
ts_decompose(USgas)Se puede apreciar en la gráfica anterior que la tendencia de la serie es bastante plana entre 2000 y 2010, y tiene un crecimiento bastante lineal. Por esa razón, la tendencia general entre 2000 y 2018 no es del todo lineal. Antes de usar la función lm, la función de regresión lineal R incorporada del paquete stats, tendremos que transformar la serie de un objeto ts a un objeto data.frame. Por lo tanto, utilizaremos la función ts_to_prophet del paquete TSstudio
USgas_df <- ts_to_prophet(USgas)La función transforma el objeto ts en dos columnas de data.frame, donde las dos columnas representan los componentes de tiempo y numéricos de la serie, respectivamente:
head(USgas_df)## ds y
## 1 2000-01-01 2510.5
## 2 2000-02-01 2330.7
## 3 2000-03-01 2050.6
## 4 2000-04-01 1783.3
## 5 2000-05-01 1632.9
## 6 2000-06-01 1513.1Ahora bien, después de transformar la serie en un objeto data.frame, se puede comenzar a crear las características de entrada de regresión. La primera característica que se define es la tendencia de la serie. Un enfoque básico para construir la variable de tendencia es indexar las observaciones de la serie en orden cronológico:
USgas_df$trend <- 1:nrow(USgas_df)La regresión de la serie con el índice de la serie proporciona una estimación del crecimiento marginal de mes a mes, ya que el índice está en orden cronológico con incrementos constantes.
La segunda característica que queremos crear es el componente estacional. Dado que queremos medir la contribución de cada unidad de frecuencia a la oscilación de la serie, utilizaremos una variable categórica para cada unidad de frecuencia. En el caso de la serie USgas, las unidades de frecuencia representan los meses del año, y, por lo tanto, crearemos una variable categórica con 12 categorías, cada categoría correspondiente a un month específico del año. Usaremos la función mes del paquete de lubridate para extraer el mes del año de la variable ds date:
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionUSgas_df$seasonal <- month(USgas_df$ds, label = T)Utilizamos la función factor para convertir la salida de la función mes en ninguna variable categórica ordenada. Revisemos ahora el marco de datos después de agregar las nuevas características:
head(USgas_df)## ds y trend seasonal
## 1 2000-01-01 2510.5 1 ene
## 2 2000-02-01 2330.7 2 feb
## 3 2000-03-01 2050.6 3 mar
## 4 2000-04-01 1783.3 4 abr
## 5 2000-05-01 1632.9 5 may
## 6 2000-06-01 1513.1 6 junPor último, antes de comenzar a retroceder la serie con esas características, dividiremos la serie en una partición de entrenamiento y prueba. Estableceremos los últimos 12 meses de la serie como una partición de prueba:
h <- 12 # setting a testing partition length
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]Ahora, al crear los marcos de datos de entrenamiento y prueba, revisemos cómo el modelo de regresión captura cada uno de los componentes por separado y todos juntos.
Modelado de la tendencia de la serie y los componentes estacionales
Primero modelaremos la tendencia de la serie haciendo retroceder la serie con la variable de tendencia, en la partición de entrenamiento:
md_trend <- lm(y ~ trend, data = train)
#Se usó la función summary para ver los detalles del modelo
summary(md_trend)##
## Call:
## lm(formula = y ~ trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -547.2 -307.4 -153.2 333.1 1052.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1751.0074 52.6435 33.26 < 2e-16 ***
## trend 2.4489 0.4021 6.09 4.86e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 394.4 on 224 degrees of freedom
## Multiple R-squared: 0.1421, Adjusted R-squared: 0.1382
## F-statistic: 37.09 on 1 and 224 DF, p-value: 4.861e-09Como puede ver en la salida de regresión anterior, el coeficiente de la variable de tendencia es estadísticamente significativo a un nivel de 0.001. Sin embargo, el R-cuadrado ajustado de la regresión es bastante bajo, lo que generalmente tiene sentido, ya que la mayor parte de la variación de la serie está relacionada con el patrón estacional como se pudo ver en las gráficas anteriores.
Para poner en contexto los datos se usará el modelo que hemos creado para predecir los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La función predict del paquete de estadísticas, como su nombre lo indica, predice los valores de un dato de entrada basado en un modelo dado.
Esto se utilizó para predecir tanto los valores ajustados como los previstos del modelo de tendencia que entrenamos antes y lo graficamos usando el paquete plotly
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)Ahora se crea la función para gráficar el modelo de salida.
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = ""),
yaxis = list(title = "Billion Cubic Feet"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Vamos a establecer las entradas de la función plot_lm con la salida del modelo:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")En general, el modelo fue capaz de capturar el movimiento general de la tendencia, sin embargo, una tendencia lineal puede no capturar la ruptura estructural de la tendencia que ocurrió alrededor de 2010.
Para el análisis de comparación, queremos medir la tasa de error del modelo tanto en el entrenamiento como en los conjuntos de prueba:
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.1644088 0.1299951El proceso de modelización y previsión del componente estacional sigue el mismo proceso que aplicamos con la tendencia, haciendo retroceder la serie con la variable estacional que creamos antes:
md_seasonal <- lm(y ~ seasonal, data = train)
#Ahora se revisaron los detalles del modelo según el siguiente código
summary(md_seasonal)##
## Call:
## lm(formula = y ~ seasonal, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -608.98 -162.34 -50.77 148.40 566.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2030.88 14.43 140.747 < 2e-16 ***
## seasonal.L -480.00 50.24 -9.554 < 2e-16 ***
## seasonal.Q 1103.83 50.17 22.000 < 2e-16 ***
## seasonal.C 72.42 50.05 1.447 0.149389
## seasonal^4 174.60 50.07 3.487 0.000592 ***
## seasonal^5 288.01 50.13 5.745 3.13e-08 ***
## seasonal^6 -44.67 50.09 -0.892 0.373460
## seasonal^7 -187.91 49.96 -3.762 0.000218 ***
## seasonal^8 84.95 49.84 1.705 0.089706 .
## seasonal^9 46.16 49.78 0.927 0.354828
## seasonal^10 77.55 49.76 1.559 0.120587
## seasonal^11 -11.09 49.75 -0.223 0.823856
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 216.9 on 214 degrees of freedom
## Multiple R-squared: 0.7521, Adjusted R-squared: 0.7394
## F-statistic: 59.04 on 11 and 214 DF, p-value: < 2.2e-16Debido a que se ha hecho retroceder la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías, que son las incrustadas con los valores de pendiente. Como se puede ver en el resumen de regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. Además, puede notar que el R-cuadrado ajustado del modelo estacional es algo más alto con respecto al modelo de tendencia (0.78 en lugar de 0.1)
Antes de trazar el modelo ajustado y pronosticar los valores con la función plot_lm, actualizaremos los valores de yhat con la función predict:
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")La gráfica anterior muestra que el modelo está haciendo un trabajo bastante bueno al capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que falta la tendencia de la serie. Antes de agregar tanto la tendencia como los componentes estacionales, calificaremos el rendimiento del modelo:
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.08628973 0.21502100La alta tasa de error en el conjunto de pruebas está relacionada con el componente de tendencia que no se incluyó en el modelo. El siguiente paso es unir los dos componentes en un modelo y pronosticar los valores de características de la serie:
md1 <- lm(y ~ seasonal + trend, data = train)
#Se revisa el resumen del contenido de la tendencia de las componentes estacionales
summary(md1)##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -514.73 -77.17 -17.70 85.80 336.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1733.7153 17.0794 101.509 < 2e-16 ***
## seasonal.L -498.1709 29.6354 -16.810 < 2e-16 ***
## seasonal.Q 1115.2951 29.5872 37.695 < 2e-16 ***
## seasonal.C 78.9835 29.5103 2.676 0.00802 **
## seasonal^4 175.6505 29.5196 5.950 1.09e-08 ***
## seasonal^5 285.0192 29.5568 9.643 < 2e-16 ***
## seasonal^6 -49.3611 29.5319 -1.671 0.09610 .
## seasonal^7 -192.3050 29.4540 -6.529 4.77e-10 ***
## seasonal^8 81.8787 29.3835 2.787 0.00581 **
## seasonal^9 44.4849 29.3480 1.516 0.13106
## seasonal^10 76.8636 29.3372 2.620 0.00943 **
## seasonal^11 -11.2755 29.3353 -0.384 0.70109
## trend 2.6182 0.1305 20.065 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared: 0.9142, Adjusted R-squared: 0.9094
## F-statistic: 189.2 on 12 and 213 DF, p-value: < 2.2e-16La regresión de la serie con los componentes de tendencia y estacional juntos proporciona una elevación adicional al R-cuadrado ajustado del modelo de 0.78 a 0.91. Esto se puede ver en la gráfica de la salida del modelo:
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Se mide acontinuación la puntuación MAPE del modelo tanto en las particiones de entrenamiento como en las de prueba:
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.04774945 0.09143290La regresión de la serie con los componentes de tendencia y estacionales proporciona un aumento significativo tanto en la calidad de ajuste del modelo como en la precisión del modelo. Sin embargo, al observar la gráfica del ajuste y pronóstico del modelo, puede notar que la tendencia del modelo es demasiado lineal y falta la ruptura estructural de la tendencia de la serie. Este es el punto en el que agregar un componente polinómico para el modelo podría proporcionar una mejora adicional para la precisión del modelo.
Una técnica simple para capturar una tendencia no lineal es agregar un componente polinómico a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Utilizaremos el argumento I, que nos permite aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Por lo tanto, usaremos este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)
#Se observa el resumen de los datos
summary(md2)##
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.882e+03 2.199e+01 85.568 < 2e-16 ***
## seasonal.L -4.917e+02 2.530e+01 -19.438 < 2e-16 ***
## seasonal.Q 1.121e+03 2.525e+01 44.381 < 2e-16 ***
## seasonal.C 8.247e+01 2.518e+01 3.275 0.00123 **
## seasonal^4 1.763e+02 2.519e+01 6.999 3.35e-11 ***
## seasonal^5 2.835e+02 2.522e+01 11.243 < 2e-16 ***
## seasonal^6 -5.175e+01 2.520e+01 -2.054 0.04123 *
## seasonal^7 -1.946e+02 2.513e+01 -7.741 3.97e-13 ***
## seasonal^8 8.030e+01 2.507e+01 3.203 0.00157 **
## seasonal^9 4.362e+01 2.504e+01 1.742 0.08293 .
## seasonal^10 7.651e+01 2.503e+01 3.057 0.00253 **
## seasonal^11 -1.137e+01 2.503e+01 -0.454 0.65005
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16La adición del polinomio de segundo grado al modelo de regresión no condujo a una mejora significativa de la bondad de ajuste del modelo. En el otro modelo, como se puede ver en la siguiente gráfica de salida del modelo, este simple cambio en la estructura del modelo nos permite capturar la ruptura estructural de la tendencia a lo largo del tiempo:
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")En el modelo que sigue la puntuación MAPE, la precisión del modelo mejoró significativamente al agregar la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas se redujo del 9,2% al 4,5%:
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2## [1] 0.03688770 0.04212618La función tslm
Hasta ahora, hemos visto el proceso manual de transformar un objeto ts en un formato de modelo de pronóstico de regresión lineal. La función tslm del paquete de pronóstico proporciona una función integrada para transformar un objeto ts en un modelo de pronóstico de regresión lineal. Con la función tslm, puede establecer el componente de regresión junto con otras características.
Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie USgas con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. Primero, dividamos la serie en particiones de entrenamiento y prueba usando la función ts_split:
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testA continuación, aplicaremos la misma fórmula que utilizamos para crear el modelo de pronóstico md2 anterior utilizando la función tslm:
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoomd3 <- tslm(train.ts ~ season + trend + I(trend^2))
summary(md3)##
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## season2 -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## season3 -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## season4 -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## season5 -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## season6 -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## season7 -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## season8 -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## season9 -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## season10 -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## season11 -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## season12 -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16Como se puede observar en la salida anterior, ambos modelos (md2 y md3) son idénticos. Hay varias ventajas de usar la función tslm, en lugar de establecer manualmente un modelo de regresión para la serie con la función lm:
- Eficiencia: no requiere transformar la serie en un objeto data.frame e ingeniería de características
- El output de salida admite toda la funcionalidad de la forecast (como las funciones de accuracy y checkresiduals) y los paquetes TSstudio (como las funciones test_forecast y plot_forecast)
Modelado de eventos individuales y eventos no estacionales
En algunos casos, los datos de series temporales pueden contener patrones inusuales que están reapareciendo con el tiempo o no. Los siguientes son ejemplos de tales eventos:
- Valores atípicos: Un solo evento o eventos que están fuera de los patrones normales de la serie.
- Ruptura estructural: Un evento significativo que cambia los patrones históricos de la serie. Un ejemplo común es un cambio en el crecimiento de la serie.
- Eventos recurrentes no estacionales: Un evento que se repite de un ciclo a otro, pero el momento en que ocurren cambia de un ciclo a otro. Un ejemplo común de tal evento son las vacaciones de Pascua, que ocurren cada año alrededor de marzo / abril.
Al no expresarse en el modelo de regresión, este tipo de eventos sesgará los coeficientes estimados, ya que el modelo ponderará esos tipos de eventos junto con los eventos regulares de la serie. El uso de variables de codificación en caliente, binarias o de indicador podría ayudar al modelo a ignorar este tipo de eventos o ajustar los coeficientes del modelo en consecuencia.
Por ejemplo, se puede observar en la gráfica de descomposición de la serie USgas mostrada anteriormente que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Si bien el crecimiento antes del año 2010 fue relativamente plano, la pendiente de la tendencia cambió después, con un crecimiento positivo. En este caso, podemos usar una variable binaria que sea igual a cero para las observaciones antes del año 2010 y un año después.
La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. En el ejemplo siguiente se muestra el proceso de creación de una variable binaria externa que es igual a 0 antes del año 2010 y a 1 posterior, utilizando la tabla USgas_df:
r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse((USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1
#Se usa la siguiiente linea de código para remodelar la serie USgas
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)
#A continuación se usa la función summary para el modelo output
summary(md3)##
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break,
## data = USgas_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -473.54 -55.44 -0.79 63.31 285.69
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.648e+03 3.166e+01 83.653 < 2e-16 ***
## season2 -3.058e+02 3.476e+01 -8.799 3.77e-16 ***
## season3 -4.857e+02 3.476e+01 -13.974 < 2e-16 ***
## season4 -9.283e+02 3.476e+01 -26.709 < 2e-16 ***
## season5 -1.109e+03 3.476e+01 -31.912 < 2e-16 ***
## season6 -1.129e+03 3.476e+01 -32.469 < 2e-16 ***
## season7 -9.591e+02 3.476e+01 -27.591 < 2e-16 ***
## season8 -9.366e+02 3.476e+01 -26.943 < 2e-16 ***
## season9 -1.141e+03 3.477e+01 -32.829 < 2e-16 ***
## season10 -1.044e+03 3.477e+01 -30.017 < 2e-16 ***
## season11 -7.927e+02 3.522e+01 -22.507 < 2e-16 ***
## season12 -2.683e+02 3.522e+01 -7.618 7.18e-13 ***
## trend -1.560e+00 4.168e-01 -3.743 0.000231 ***
## I(trend^2) 1.869e-02 1.689e-03 11.063 < 2e-16 ***
## s_break NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.9 on 224 degrees of freedom
## Multiple R-squared: 0.9412, Adjusted R-squared: 0.9377
## F-statistic: 275.6 on 13 and 224 DF, p-value: < 2.2e-16Como se puede ver en el resumen del modelo anterior, la variable de rotura estructural es estadísticamente significativa, con un nivel de 0,01. Del mismo modo, en el caso de valores atípicos o días festivos, la codificación en caliente se puede aplicar estableciendo una variable binaria que sea igual a 1 cada vez que ocurra un evento atípico o no estacional, y 0 en caso contrario.
Previsión de una serie con componentes multiestacionales: un caso de estudio
Una de las principales ventajas del modelo de regresión, frente a los modelos tradicionales de series temporales como ARIMA o Holt-Winters, es que proporciona una amplia gama de opciones de personalización y nos permite modelar y pronosticar datos de series temporales complejas como series con multiestacionalidad.
En los siguientes ejemplos, se usó la serie UKgrid para demostrar el enfoque de pronóstico de una serie multiestacional con un modelo de regresión lineal
La serie UKgrid
La serie UKgrid representa la demanda nacional de electricidad en el Reino Unido, y está disponible en el paquete UKgrid. Esta serie representa una serie temporal de datos de alta frecuencia con frecuencia de media hora. Utilizaremos la función extract_grid del paquete UKgrid para definir la serie, las características principales (por ejemplo, formato de datos, variables, frecuencia, etc.). Esta función de transformación nos permite agregar la frecuencia de la serie de media hora a una frecuencia más baja, como horaria, diaria o mensual. Como nuestro objetivo aquí es pronosticar la demanda diaria en los próximos 365 días, estableceremos la serie a la frecuencia diaria utilizando la estructura data.frame
library(UKgrid)
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
head(UKdaily)## TIMESTAMP ND
## 1 2005-04-01 1920069
## 2 2005-04-02 1674699
## 3 2005-04-03 1631352
## 4 2005-04-04 1916693
## 5 2005-04-05 1952082
## 6 2005-04-06 1964584ts_plot(UKdaily, title = "The UK National Demand for Electricity", Ytitle = "MW", Xtitle = "Year")- Diario: Un ciclo de 365 días al año
- Día de la semana: Un ciclo de 7 días
- Mensual: Afectado por el clima
La evidencia de esos patrones se puede ver en el siguiente mapa de calor de la serie desde 2016 utilizando la función ts_heatmap del paquete TSstudio:
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")Preprocesamiento e ingeniería de características de la serie UKdaily
Para capturar los componentes estacionales de la serie, estableceremos la serie como frecuencia diaria y crearemos las siguientes dos características:
-Indicador de día de la semana -Indicador del mes del año
Además, como es razonable suponer (confirmaremos esta suposición con la función ACF una vez que hayamos transformado la serie en un objeto ts) que la serie tiene una fuerte correlación con los retrasos estacionales, crearemos una variable de retraso con un retraso de 365 observaciones. Utilizaremos el paquete dplyr para crear esas características:
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionUKdaily <- UKdaily %>%
mutate(wday = wday(TIMESTAMP, label = TRUE),
month = month(TIMESTAMP, label = TRUE),
lag365 = dplyr::lag(ND, 365)) %>%
filter(!is.na(lag365)) %>%
arrange(TIMESTAMP)
#
str(UKdaily)## 'data.frame': 4939 obs. of 5 variables:
## $ TIMESTAMP: Date, format: "2006-04-01" "2006-04-02" ...
## $ ND : int 1718405 1691341 1960325 2023886 2026204 2008422 1981175 1770440 1749715 2012865 ...
## $ wday : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 7 1 2 3 4 5 6 7 1 2 ...
## $ month : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 4 4 4 4 4 4 4 4 4 4 ...
## $ lag365 : int 1920069 1674699 1631352 1916693 1952082 1964584 1990895 2003982 1811436 1684720 ...Como la entrada de la función tslm debe estar en formato ts (al menos para la serie), convertiremos la serie a un objeto ts. Usaremos la primera marca de tiempo de la serie y las funciones year and yday (el día del año) del paquete lubridate para establecer el punto de partida del objeto:
start_date <- min(UKdaily$TIMESTAMP)
UK_ts <- ts(UKdaily$ND,
start = c(year(start_date), yday(start_date)),
frequency = 365)Después de transformar la serie en un objeto ts, podemos volver atrás y confirmar la suposición que hicimos sobre el nivel de correlación entre la serie y sus retrasos estacionales con la función ts_acf (una versión interactiva de la función acf del paquete TSstudio). Revisaremos la correlación de la serie con sus rezagos de los últimos cuatro años.
acf(UK_ts, lag.max = 365 * 4)
Ahora, después de haber creado las nuevas características para el modelo y establecido el objeto ts, estamos listos para dividir la serie de entrada y el objeto de características externas correspondiente que creamos (UKdaily) en una partición de entrenamiento y prueba. Como nuestro objetivo es pronosticar las próximas 365 observaciones, y la longitud de la serie es lo suficientemente grande (2.540 observaciones), podemos permitirnos establecer la partición de prueba a la longitud del horizonte de pronóstico: 365 observaciones. Estableceremos h, una variable indicadora, en 365 y la usaremos para definir las particiones y, más adelante, el horizonte de pronóstico:
h <- 365
#Como antes, dividiremos la serie en una partición de entrenamiento y prueba con la función ts_split:
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$testDe manera similar, tenemos que dividir las características que creamos para el modelo de regresión (las características estacionales y de retraso) en una partición de entrenamiento y prueba siguiendo exactamente el mismo orden que usamos para el objeto ts correspondiente. Usaremos la funcionalidad de índice data.frame para establecer la tabla UKdaily en particiones de entrenamiento y prueba:
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]Entrenamiento y prueba del modelo de pronóstico
Después de haber creado las diferentes características para el modelo, estamos listos para comenzar el proceso de entrenamiento del modelo de pronóstico. Utilizaremos la partición de entrenamiento y comenzaremos a entrenar los siguientes tres modelos:
- Modelo de línea base: Regresión de la serie con el componente estacional y de tendencia utilizando las características integradas de la función tslm. Como establecemos la frecuencia de la serie en 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
- Modelo multiestacional: Sumando los indicadores de día de la semana y mes del año para capturar la multiestacionalidad de la serie.
- Un modelo multiestacional con un desfase estacional: Utilizando, además de los indicadores estacionales, la variable de rezago estacional.
La comparación de estos tres modelos se basará en los siguientes criterios: - Rendimiento del modelo en el conjunto de entrenamiento y pruebas utilizando la puntuación MAPE - Visualización de los valores ajustados y pronosticados frente a los valores reales de la serie mediante la función test_forecast
Comenzaremos con el modelo de referencia, haciendo retroceder la serie con sus componentes estacionales y de tendencia:
md_tslm1 <- tslm(train_ts ~ season + trend)A continuación, utilizaremos el modelo entrenado, md_tslm1, para pronosticar los próximos 365 días de la serie, correspondientes a las observaciones de la partición de prueba, utilizando la función de pronóstico:
fc_tslm1 <- forecast(md_tslm1, h = h)Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y prueba utilizando la función test_forecast:
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)Podemos observar en la gráfica de rendimiento anterior que el modelo de referencia está haciendo un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año. Por otro lado, no logra captar la oscilación que se relaciona con el día de la semana. La puntuación MAPE del modelo, como podemos ver en la salida de la siguiente función de precisión, es de 6,09% y 7,77% en las particiones de entrenamiento y prueba, respectivamente:
accuracy(fc_tslm1, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.030002e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set -1.740215e+04 123156.6 96785.27 -1.8735336 7.160573 0.8112374
## ACF1 Theil's U
## Training set 0.5277884 NA
## Test set 0.5106681 1.071899Ahora intentaremos mejorar la precisión del modelo, agregando las características del día de la semana y el mes del año al modelo:
md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)Como ahora estamos utilizando características de una fuente de datos externa, tenemos que especificar los datos de entrada con el argumento de datos. Ahora repetiremos el mismo proceso que antes, utilizando un pronóstico con el modelo entrenado y evaluando el rendimiento del modelo:
fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm2,
test = test_ts)Mirando la gráfica de rendimiento anterior del segundo modelo, podemos observar la contribución de las características estacionales en el pronóstico, ya que el segundo modelo pudo capturar tanto la tendencia como los patrones multiestacionales de la serie. Esto también se puede observar en el puntaje MAPE del modelo, que cayó a 2.87% y 5.23% en las particiones de entrenamiento y prueba, respectivamente:
accuracy(fc_tslm2, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.812170e-13 70245.98 52146.79 -0.1738708 3.167605 0.4370853
## Test set -1.764563e+04 80711.71 65373.21 -1.3715505 4.682071 0.5479470
## ACF1 Theil's U
## Training set 0.7513664 NA
## Test set 0.6075598 0.68445Por último, pero no menos importante, se agrera la variable lag al modelo y se repite el mismo proceso que antes:
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365, data = train_df)
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)Es difícil observar a partir de la gráfica de rendimiento del tercer modelo, si hay una diferencia significativa con respecto al segundo modelo. Por lo tanto, revisemos el puntaje MAPE del tercer modelo:
accuracy(fc_tslm3, test_ts) ## ME RMSE MAE MPE MAPE MASE
## Training set -1.175109e-11 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
## ACF1 Theil's U
## Training set 0.7500146 NA
## Test set 0.6094666 0.6925767Los resultados del tercer modelo muestran una pequeña mejora en la precisión del modelo con un 2,81% en el conjunto de entrenamiento y un 5,07% en el conjunto de pruebas.
Selección de modelos
Ahora es el momento de seleccionar un modelo. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Como tanto el segundo como el tercer modelo lograron una puntuación MAPE bastante similar, con una pequeña ventaja para el tercer modelo, deberíamos preguntarnos si una mejora de menos del 0,2% en la tasa de error del conjunto de pruebas vale la pena el costo de usar la variable lag (es decir, la pérdida de 365 observaciones y el costo adicional de un grado de libertad para el modelo).
No hay una respuesta directa a esta pregunta. Además, depende del objetivo del pronóstico. Se recomienda que considere la siguiente prueba:
La primera pregunta que se debe hacer en este caso: ¿Es la variable de retraso estadísticamente significativa? Si la variable no es estadísticamente significativa, no tiene sentido continuar la discusión, y es mejor abandonar la variable. En el caso del tercer modelo, podemos utilizar la función de resumen para observar el nivel de significancia de la variable lag365:
summary(md_tslm3)$coefficients %>% tail(1)## Estimate Std. Error t value Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05El valor p de la variable lag365 indicó que la variable es estadísticamente significativa a un nivel de 0,001. Del mismo modo, podemos aplicar una única prueba ANOVA con la función anova del paquete stats, y comprobar si la variable adicional es significativa
anova(md_tslm3)## Analysis of Variance Table
##
## Response: train_ts
## Df Sum Sq Mean Sq F value Pr(>F)
## season 364 1.4279e+14 3.9227e+11 73.5348 < 2.2e-16 ***
## trend 1 7.2634e+13 7.2634e+13 13615.7078 < 2.2e-16 ***
## wday 6 4.5009e+13 7.5016e+12 1406.2214 < 2.2e-16 ***
## month 11 1.3721e+11 1.2473e+10 2.3382 0.007266 **
## lag365 1 8.1432e+10 8.1432e+10 15.2650 9.49e-05 ***
## Residuals 4190 2.2352e+13 5.3345e+09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Esta prueba ha arrogado que la variable estadísticamente significativa
Podría darse el caso de que el tercer modelo sea más preciso solo por casualidad y no porque la variable lag contribuya a la precisión del modelo, ya que la diferencia es relativamente pequeña. Por lo tanto, el backtesting de ambos modelos podría ayudar a validar si la contribución de la variable lag es consistente durante varios períodos de prueba. Te dejaré a ti realizar backtesting para ambos modelos como un ejercicio divertido.
En aras de la simplicidad, iremos con los criterios de precisión y seleccionaremos el tercer modelo para pronosticar la serie. El siguiente paso es volver a entrenar el modelo en todas las series y pronosticar los próximos 365 días:
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365,
data = UKdaily)Análisis de residuos
Justo antes de finalizar el pronóstico, hagamos un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals:
checkresiduals(final_md)
##
## Breusch-Godfrey test for serial correlation of order up to 730
##
## data: Residuals from Linear regression model
## LM test = 3301.1, df = 730, p-value < 2.2e-16Como puede ver en el gráfico de resumen residual anterior, los residuos no son ruido blanco, ya que existe cierta autocorrelación entre las series residuales y sus retrasos. Esto es técnicamente una indicación de que el modelo no capturó todos los patrones o información que existe en la serie. Una forma de abordar este problema es identificar variables adicionales que puedan explicar la variación en los residuos. El principal desafío con este enfoque es que es difícil identificar variables externas que puedan explicar la variación de los residuos, y también sean factibles de pronosticar. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, sin embargo, es difícil predecir esos patrones climáticos un año antes.
Un enfoque alternativo, cuando los patrones sobrantes en los residuos del modelo es tratar los residuos del modelo como datos de series temporales separados y modelarlos con otro modelo de pronóstico de series temporales.
Finalización del pronóstico
Finalicemos el proceso y utilicemos el modelo entrenado seleccionado para pronosticar las futuras observaciones 365. Dado que utilizamos variables externas con la función tslm, tendremos que generar sus valores futuros. Esto es relativamente simple, ya que utilizamos variables deterministas. Por lo tanto, crearemos un objeto data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras. Un enfoque simplista es crear primero las fechas correspondientes de las observaciones pronosticadas, y luego extraer de este objeto el día de la semana y el mes del año:
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1),
by = "day",
length.out = h))A continuación, se utilizaron la variable de fecha para crear las variables wday y month con el paquete de lubricante:
UK_fc_df$wday <- factor(lubridate::wday(UK_fc_df$date, label = TRUE), ordered = FALSE)
UK_fc_df$month <- factor(month(UK_fc_df$date, label = TRUE), ordered = FALSE)
UK_fc_df$lag365 <- tail(UKdaily$ND, h)Por último se trazó el pronóstico final con la función plot_forecast del paquete TSstudio:
#También se creó el forcast final
UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)
plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")Resumen
Gracias al estudio de este capítulo del libro aprendimos como representar las aplicaciones de pronóstico del modelo de regresión lineal. Aunque el modelo de regresión lineal no fue diseñado para manejar datos de series temporales, gracias a una ingeniería de características simple podemos transformar un problema de pronóstico en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal con respecto a otros modelos tradicionales de series temporales es la capacidad del modelo para incorporar variables y factores externos. Sin embargo, este modelo puede manejar series temporales con patrones de multiestacionalidad, como se observó con el pronóstico de demanda de electricidad del Reino Unido.