El presente trabajo es un resumen de los aspectos tratados en el texto Hands On Time Series Analysis with R en el cual se replicaron los códigos para la enseñanza de la estadística implicada, por: Danna Patiño Rueda y Daniel Rivera Bolaño.

¿Regresión lineal? ¿Qué es?

La regresión lineal es una técnica de modelado estadístico que permite crear un modelo líneal, el cual describe una relación lineal entre una variable dependiente, \(y\), y una variable(s) independiente(s), \(x\). Así, la ecuación general que corresponde a dicho modelo de regresión lineal es:

Para el caso de una sola variable independiente: \[Y_i= \beta_0+\beta_1 \cdot X_{1,i}+ \in_i \] Para n variables independientes: \[Y_i =\beta_0 +\beta_1 X_{1,i}+\beta_2X_{2,i}+...+\beta_nX_{n,i}+\in_i \] Donde:

Para estimar los coeficientes del modelo, debe realizarse lo siguiente:

El criterio de valoración más usual 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:

Estimación de coeficientes con el método OLS

Características de la ingeniería de los componentes de la serie

Para crear las entradas de regresión que representan la tendencia de la serie y los componentes estacionales, primero tenemos que entender su estructura. Es por esto que, 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 ha transformado en una estructura aditiva, la modificación a una formación de regresión lineal es sencilla y se rige por 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, las principales características de la serie se dan 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 10

Como se puede ver en la trama de la serie USgas, esta es una serie mensual con un fuerte componente estacional mensual y una línea de tendencia bastante estable. Podemos examinar la estructura de los componentes de la serie con la función ts_decompose más adelante:

ts_decompose(USgas)

De la gráfica anterior, se puede apreciar que la tendencia de la serie es bastante plana y tiene un crecimiento lineal entre 2000 y 2010. Por lo anterior, 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, se tendrá que modificar 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, respectivamente, las dos columnas representan los componentes de tiempo y numéricos de la serie:

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.1

Después de transformar la serie en un objeto data.frame, se comienzan 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, puesto que el índice está en orden cronológico con incrementos constantes.

La segunda característica es el componente estacional. Dado que se quiere medir la contribución de cada unidad de frecuencia a la oscilación de la serie, se utilizará 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, se creará una variable categórica con 12 categorías, cada categoría correspondiente a un month específico del año. Se usará 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, union
USgas_df$seasonal <- month(USgas_df$ds, label = T)

Se utilizó 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      jun

Por último, antes de comenzar a retroceder la serie con esas características, dividiremos la serie en una partición de entrenamiento y prueba. Se establecerán 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 se modelará la tendencia de la serie haciendo retrocederla 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-09

Como se 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 tiene sentido, puesto que la mayor parte de la variación de la serie está relacionada con el patrón estacional justo como se pudo ver en las gráficas anteriores.

Para poner en contexto los datos se usará el modelo que se ha 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)

Ya se puede crear 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':
## 
##     layout
plot_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)
}

Establecemos 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")

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ó aproximadamente en 2010.

Para el análisis de comparación, se quiere 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.1299951

El proceso de modelización y previsión del componente estacional se rige por 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-16

Dado 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 implementadas 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, se puede notar que el R-cuadrado ajustado del modelo estacional es más alto con respecto al modelo de tendencia (0.78 en lugar de 0.1)

Se actualizarán los valores de yhat con la función predict, antes de trazar el modelo ajustado y pronosticar los valores con la función plot_lm:

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")

En la gráfica anterior se puede observar que el modelo está haciendo un buen trabajo al capturar la estructura del patrón estacional de la serie. No obstante, se puede evidenciar que la tendencia de la serie es ausente. Y es por esto que, antes de agregar tanto la tendencia como los componentes estacionales, se calificará 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.21502100

La alta tasa de error en el conjunto de pruebas se debe a que no se incluyó el componente de tendencia en el modelo. Es por esto que ahora, se deben 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-16

La regresión de la serie con los componentes de tendencia y estacional juntos, proporcionan una elevación adicional al R-cuadrado ajustado del modelo de 0.78 a 0.91. Esto se puede observar 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 medirá 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.09143290

La 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 de este. No obstante, al observar la gráfica del ajuste y pronóstico del modelo, se puede notar que la tendencia del modelo es 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 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. Se utilizará el argumento I para aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Es por esto que se usará dicho 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-16

La 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 podrá evidenciar en la siguiente gráfica de salida del modelo, este simple cambio 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, puesto 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.04212618

Función tslm

Se ha visto el proceso manual de transformar un objeto ts en un formato de modelo de pronóstico de regresión lineal. Ahora bien, 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, se puede establecer el componente de regresión junto con otras características.

Se repetirá el ejemplo anterior y se pronosticarán 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. Para esto se dividirá 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$test

Luego, se aplicará 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 zoo
md3 <- 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-16

Como se pudo observar en la salida anterior, ambos modelos (md2 y md3) son idénticos. 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 es necesario 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 o no reapareciendo con el tiempo, por ejemplo:

Al no expresarse en el modelo de regresión, este tipo de eventos sesgará los coeficientes estimados, puesto 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.

Por ejemplo, en la gráfica de descomposición de la serie USgas mostrada anteriormente se puede ver 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ó con un crecimiento positivo. Para este caso, se puede usar una variable binaria que sea igual a cero para las observaciones antes y un año después del año 2010.

La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. En el siguiente ejemplo se muestra el proceso de creación de una variable binaria externa que es igual a 0 y posterior a 1 del año 2010, respectivamente, 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-16

Como se pudo observar en el resumen del modelo anterior, la variable de rotura estructural es estadísticamente significativa, con un nivel de 0,01. Así mismo, en el caso de valores atípicos (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, permitiendonos modelar y pronosticar datos de series temporales complejas como series con multiestacionalidad.

En los proximos ejemplos, se utilizó la serie UKgrid para demostrar el enfoque de pronóstico de una serie multiestacional con un modelo de regresión lineal.

Serie UKgrid

La serie UKgrid representa la demanda nacional de electricidad en el Reino Unido, y está disponible en el paquete UKgrid. Esta representa una serie temporal de datos de alta frecuencia cada media hora. Se utilizará la función extract_grid del paquete UKgrid para definir la serie, las características principales (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. Dado que 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 1964584
ts_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

Esos patrones se pueden evidenciar 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, se establecerá la serie como frecuencia diaria y se crearán las siguientes características:

-Indicador de día de la semana -Indicador del mes del año

Además, se supondrá (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, se creará una variable de retraso con un retraso de 365 observaciones. Se utilizará 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, union
UKdaily <- 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 (para la serie), se convertirá la serie a un objeto ts. Se usará 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)

Luego de transformar la serie a un objeto ts, se puede volver atrás y confirmar la suposición que se hizo 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). Se revisará la correlación de la serie con sus rezagos de los últimos cuatro años.

acf(UK_ts, lag.max = 365 * 4)

Después de haber creado las nuevas características para el modelo y establecido el objeto ts, se podrá 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), se puede establecer la partición de prueba a la longitud del horizonte de pronóstico: 365 observaciones. Se establecerá h (una variable indicadora) en 365 y la se usará 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$test

Se tendrá 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. Se usará 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

Se iniciará el proceso de entrenamiento del modelo de pronóstico. Se utilizará la partición de entrenamiento y se comenzará 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 se estableció 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 modelos se darpa por medio de 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

Se iniciará con el modelo de referencia, haciendo retroceder la serie con sus componentes estacionales y de tendencia:

md_tslm1 <- tslm(train_ts ~ season + trend)

Luego, se utilizará 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)

En la gráfica de rendimiento anterior se pudo observar que el modelo de referencia está realizando un buen 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 observar 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.071899

Ahora mejoraremos la precisión del modelo, agregando las características del día de la semana y el mes del año:

md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)

Dado que ahora estamos utilizando características de una fuente de datos externa, tenemos que especificar los datos de entrada con el argumento de datos. Repetiremos el mismo proceso que antes, utilizando un pronóstico con el modelo entrenado y evaluando el rendimiento de este:

fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
              forecast.obj = fc_tslm2,
              test = test_ts)

Al observar la gráfica de rendimiento anterior del segundo modelo, se puede evidenciar la contribución de las características estacionales en el pronóstico, puesto que el segundo modelo pudo aprehender 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.68445

Por último, se agrerá la variable lag al modelo y se repetirá 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 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.6925767

Los resultados del tercer modelo muestran una pequeña mejora en la precisión con un 2,81% en el conjunto de entrenamiento y un 5,07% en el conjunto de pruebas.

Selección de modelos

Ya se tiene en mente 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. Por lo que, se recomienda que considerar la siguiente prueba:

Primero, debe preguntarse: ¿Es la variable de retraso estadísticamente significativa? Si la variable no es estadísticamente significativa 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-05

El 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 ' ' 1

Esta prueba ha arrojado que la variable es 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 considerablemente 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.

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

Antes de finalizar el pronóstico, se hará 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-16

Como se puede observar 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 una indicación de que el modelo no capturó todos los patrones o información que existen en la serie. Para abordar este problema se pueden 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 que sean factibles de pronosticar. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, no obstante, es difícil predecir esos patrones climáticos un año antes.

Otro enfoque es cuando los patrones sobrantes en los residuos del modelo tratan dichos residuos como datos de series temporales separados y modelarlos con otro modelo de pronóstico de series temporales.

Finalización del pronóstico

Para finalizar, utilicemos el modelo entrenado seleccionado para pronosticar las futuras observaciones 365. Puesto que utilizamos variables externas con la función tslm, tendremos que generar sus valores futuros. Para esto utilizaremos 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. Se pueden crear primero las fechas correspondientes de las observaciones pronosticadas, y luego extraer 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))

Luego, 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 y, dada una ingeniería de características simple, podemos transformar un problema de pronóstico en un problema de regresión lineal. Una gran ventaja de este modelo (regresión lineal) con respecto a otros modelos tradicionales de series temporales es la capacidad del modelo para incorporar variables y factores externos. No obstante, como se observó con el pronóstico de demanda de electricidad del Reino Unido, este modelo puede manejar series temporales con patrones de multiestacionalidad.