El siguiente trabajo se basa en un resumen del capitulo “Pronóstico de Regresión Lineal” (Cap. 9) del libro “Hands On Time Series Analysis with R” donde tomamos como base estos códigos para la realización del informe, presentado por Michael Albornoz y David Cotes.

¿Qué es la regresión lineal?

El modelo de regresión lineal busca identificar y cuantificar la relación entre \(y\) (variable dependiente) y \(x\). (independiente) sea univariante o multivariante. Este modelo se puede aplicar desde la inferencia causal, el análisis predictivo hasta la previsión de series temporales. Generalmente a \(y\) se le conoce como variable de respuesta y a \(x\) como variable predictora, impulsora o regresora. Siguiendo la misma idea, este modelo se expresaría como una combinación lineal de las variables independientes. En esta aplicación nos enfocaremos en los métodos y enfoques para pronosticar datos de series de tiempo con regresión lineal, descomponiendo y realizando un pronóstico de los componentes de la serie (la tendencia y los patrones estacionales), en manejo de eventos especiales (valores atípicos y días festivos) y el uso de variables externas como regresores.

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:

La estimación de los coeficientes del modelo se basa en los siguientes dos pasos:

El enfoque de estimación más común es aplicar el Mínimos cuadrados ordinarios (MCO) para minimizar los residuos de la suma de cuadrados 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

Antes de la creación de entradas de regresión para representar la tendencia de la serie y los componentes estacionales, debemos comprender su estructura. Para los siguientes ejemplos, se demuestran los procesos de creación de nuevas funciones a partir de series existentes implementando la serie 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 10

Como podemos observar 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. Para explorar la estructura de los componentes de la serie, usamos la función ts_decompose:

ts_decompose(USgas)

Observando la gráfica anterior, la tendencia de la serie es bastante plana entre 2000 y 2010, y creciendo linealmente en el futuro. Dicho esto, la tendencia general entre 2000 y 2018 no es estrictamente lineal. Antes de utilizar lm, la función de regresión lineal R incorporada en el paquete stats, tendremos que transformar la serie de un objetor ts a un objeto data frame. Por lo tanto, utilizaremos 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.1

Seguidamente, luego de transformar la serie en un objeto data frame, comenzamos creando las características de entrada de regresión. Primeramente, definiendo la primera característica de la serie, tenemos un enfoque básico para construir la variable de tendencia indexando 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, union
USgas_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      Jan
## 2 2000-02-01 2330.7     2      Feb
## 3 2000-03-01 2050.6     3      Mar
## 4 2000-04-01 1783.3     4      Apr
## 5 2000-05-01 1632.9     5      May
## 6 2000-06-01 1513.1     6      Jun

Finalmente, para retroceder la serie con dichas características, dividimos la serie en una partición de entrenamiento y prueba, estableciendo los últimos 12 meses de la serie.

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.

Modelización 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-09

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

Poniendo en contexto los datos, se usa el modelo creado para predecir los valores ajustados en la petición de entrenamiento y los valores pronosticados en la partición de prueba. La función predict del paquete de stats, utiliza un modelo como base para predecir valores de un dato de entrada.

Esto con el fin de predecir no solo los valores ajustados sino también los previstos del modelo de tendencia que entrenamos antes y lo graficamos usando el paquete plothy.

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':
## 
##     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)
}

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

Generalmente, 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 como anteriormente se dio 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.1299951

El proceso de modelización y previsión del componente estacional sifue estando en 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

Debido 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 excelente capturando la estructura del patrón estacional de la serie. Sin embargo, podemos observar que falta la tendencia de la serie. Antes de agregar tanto la tendencia como los componentes estacionales, se debe calificar la eficiencia 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 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-16

La 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.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 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-16

La adición del polinomio de segundo grado al modelo de regresión no condujo a una mejoría de la bondad de ajuste del modelo. En el otro modelo, notando 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, sí hubo una mejoría significativa al agregar la tendencia polinómica al modelo de regresión, esto se dio gracias a que el error en el conjunto 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

La función TSLM

En lo que va de esta aplicación, hemos experimentado sobre 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. Al utilizar la función ts, esta configurará el componente de regesión junto con otras funciones.

Ahora se repetirá el ejemplo anterior y pronosticaremos las últimas 12 observaciones de USgas serie con el tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. Primero, dividamos la serie en particiones de entrenamiento y prueba usando elts_split función:

USgas_split <- ts_split(USgas, sample.out = h)

train.ts <- USgas_split$train

test.ts <- USgas_split$test

A 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 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 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:

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-16

Como se puede ver 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 o días festivos, la codificación se puede aplicar estableciendo una variable binaria que sea 0 en el caso contrario o igual a 1 cada vez que ocurra un evento atípico o no estacional.

Pronóstico de una serie con componentes de multiestacionalidad: un estudio de caso

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

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

Ahora, se creará una variable de retraso de 365 observaciones utilizando el paquete dplyr, ya que la serie está relacionada a retrasos estacionales. Todo esto para confirmar la función ACF cuando se haya transformado la serie a un objeto ts.

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 "Sun"<"Mon"<"Tue"<..: 7 1 2 3 4 5 6 7 1 2 ...
##  $ month    : Ord.factor w/ 12 levels "Jan"<"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)

Luego de transformar la serie en un objeto ts, se puede rebobinar 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_act (versión interactiva de la función act del paquete TSstudio). Se revisará la correlación de la serie con sus rezagos de los últimos 4 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, está todo listo para dividir la serie de entrada y el objeto de características externas correspondiente que creamos (Reino Unido diario) en una partición de entrenamiento y prueba. Ya que el objetivo es pronosticar las próximas 365 observaciones y la longitud de la serie es lo suficientemente grande (2540 observaciones), podemos darnos el lujo de establecer la partición de prueba en la longitud del horizonte de pronóstico: 365 observaciones. Estableceremos una variable indicadora, para 365 y utilizarlo para definir las particiones y, posteriormente, el horizonte de previsión:

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

De 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

Luego de crear las diferentes funciones para el modelo, se comienza el proceso de capacitación del modelo de pronóstico. Se usa la partición de entrenamiento y comenzaremos a entrenar los siguientes tres modelos:

  • Modelo de referencia: Regresión de la serie con el componente estacional y de tendencia utilizando las características integradas de tslm. Conforme fijamos la frecuencia de la serie en 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
  • Modelo multiestacional: Se suman indicadores de día de la semana y mes del año para captar 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)

Nótese en el gráfico 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 capturar la oscilación relacionada con el día de la semana. La puntuación MAPE del modelo es 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 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)

En el gráfico de rendimiento anterior del segundo modelo, se observa 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 la puntuación del modelo MAPE, que se redujo 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, 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.6925767

Los 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

Este es el momento para seleccionar un modelo, es evidente que el segundo y tercer modelo funcionan mejor que el primero, así como el segundo y tercer modelo logran una puntuación MAPE muy similar, dejando al tercer modelo con una pequeña ventaja hay que preguntarse si una mejora de menos del 0,2 % en la tasa de error del conjunto de prueba vale la pena el costo de usar el variable de rezago (es decir, la pérdida de 365 observaciones y el costo adicional de un grado de libertad para el modelo).

La respuesta a esto depende, no hay una respuesta directa. Depende del objetivo del pronóstico, por lo que recomienda considerar la siguiente prueba:

En primer lugar, hay que preguntarse: ¿Es la variable de rezago estadísticamente significativa? Si la variable no es estadísticamente significativa, no tiene sentido continuar la discusión y es mejor descartar la variable. En el caso del tercer modelo, puede utilizarse la función de resumen para observar el nivel de significación 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

La prueba ANOVA también indicó que la variable lag365 es estadísticamente significativa

Prueba retrospectiva: Puede darse el caso de que el tercer modelo sea más preciso sólo por casualidad y no porque la variable de retraso contribuya a la precisión del modelo, ya que la diferencia es relativamente pequeña. El backtesting de ambos modelos podría ayudar a validar si la contribución de la variable de rezago es consistente durante varios períodos de prueba.

En aras de la simplicidad, se seguirá con los criterios de precisión y seleccionar 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-16

Como se observa en el gráfico del resumen residual, estos no son de ruido blanco, ya que existe cierta autocorrelación entre las series de residuos y sus retrasos. Significa que técnicamente hay una indicación de que el modelo no capturó todos los patrones o la información que existe en la serie. Una forma de abordar este inconveniente es identificar variables adicionales que puedan explicar la variación en los residuos. El primer reto con este enfoque es la dificultad de identificar variables externas que puedan explicar la variación de los residuos y que también sean factibles de pronosticar. Un ejemplo de esto es suponer que los patrones climáticos afectan la demanda de electricidad, pero es complicado 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

Finalizando el proceso se usa el modelo entrenado seleccionado para pronosticar las 365 observaciones futuras. Al usar variables externas con la función tslm, se tendrá que generar sus valores futuros. Esto es relativamente simple, ya que usamos variables deterministas. Así que se creará un objeto data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras.

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

En el estudio de este capítulo se presentaron las aplicaciones de pronóstico del modelo de regresión lineal. Aunque este modelo no fue diseñado para manejar datos de series de tiempo, puede transformarse un problema de pronóstico en un problema de regresión lineal. Nosotros como grupo aprendimos las ventajas que esto conlleva, como identificar que el modelo de regresión lineal tiene la capacidad de incorporar variables y factores externos a diferencia de los modelos tradicionales, y aunque este modelo solo puede manejar series temporales con patrones de multiestacionalidad, como se observó con el pronóstico de demanda de electricidad presentado anteriormente.