#En este informe se presentara un resumen del texto Hands On Time Series Analysis with R de manera detallada en el cual se aplicaron los codigos que se aprendieron.

“El modelo de regresión lineal es uno de los modelos más frecuentemente usados para llevar a cabo la identificación y cuantificación de la relación que hay entre una variable dependiente y una variable independiente, tanto única como múltiple”

¿Qué es la regresión lineal?

“El uso primordial y principal del modelo de regresión lineal es cuantificar la conexión entre la variable dependiente Y (también llamada variable de respuesta), y la variable independiente X (también llamada propulsor) de manera lineal. Una relación lineal entre las variables dependiente e independiente se pueden dividir en las siguientes ecuaciones:”

“Cuando solo hay una variable independiente:”

\[Y_i= \beta_0+\beta_1 \cdot X_{1,i}+ \in_i \] “Cuando hay un número de N variables independientes:” \[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 dos pasos siguientes:

El enfoque de estimación más frecuente es aplicar los mínimos cuadrados ordinarios (Ordinary Least Squares) (OLS) como una técnica de optimización para minimizar la suma de los cuadrados de los residuos (\(\sum_{i=1}^{N}\in_i^2\)).

El cuadrado de los residuos conlleva dos efectos:

Coefficients estimation with the OLS method (Estimación de coeficientes con el método OLS)

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

El OLS es un método de optimización simple que se basa en algebra y cálculo lineal básico. El objetivo de OLS es distinguir los coeficientes que minimizan la suma de cuadrados de los residuos.

Antes de crear las entradas de regresión que representan la tendencia de la serie y los componentes estacionales, es fundamental comprender su estructura. Por ende, el libro guía demuestra el proceso de creación de nuevas funciones a partir de series existentes, tomadas de una base de datos llamada USgas.

Se grafica la serie con la función ts_plot, en base a las series de TSstudio:

library(TSstudio)

data("USgas")

ts_plot(USgas,
        title = "US Monthly Natural Gas consumption",
        Ytitle = "Billion Cubic Feet",
        Xtitle = "Year")

Presenta 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 se puede evidenciar en la gráfica de la serie, USgas es una serie mensual con un fuerte componente estacional mensual, y una línea de tendencia bastante estable. Los componentes de la serie pueden ser más explorados con la función ts_decompose:

ts_decompose(USgas)

En la anterior gráfica, se evidencia como la tendencia de la serie es altamente plana entre 2000 y 2010, además de tener un crecimiento lineal en el futuro. Por ende, la tendencia general entre el 2000 y 2018 no es completamente lineal.

Antes de utilizar la función lm, la función de regresión lineal R introducida del paquete de stats, se debe transformar la serie de un objeto ts a un objeto data.frame. Por lo tanto, se utiliza la función ts_to_prophet del paquete de TSstudio:

USgas_df <- ts_to_prophet(USgas)

La función transforma el objeto ts en donde columnas de data.frame, donde estas dos columnas representan el tiempo y los componentes 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

Posteriormente a transformar la serie en un objeto de marco de datos, es posible empezar a crear las características de entrada de la regresión.

La primera característica que se elabora es la tendencia de la serie. Un enfoque básico para llevar a cabo la construcción de la variable de tendencia es indexar el orden cronológico de las observaciones de la serie:

USgas_df$trend <- 1:nrow(USgas_df)

La regresión de la serie con el índice de serie proporciona una estimación del crecimiento marginal que hay mes a mes, debido a que el índice está en orden cronológico con incrementos constantes.

La segunda característica que se requiere elaborar es el componente estacional. Debido a que se requiere 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. Para el caso de la serie USgas, las unidades de frecuencia indican los meses del año, por ende, se crea una variable categórica con 12 categorías, donde cada categoría representa un mes específico del año. Se usa la función mes del paquete lubridate para poder 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 usa la función de factor para convertir la salida de la función mes en una variable categórica no ordenada. Posterior a agregar las nuevas funciones, se revisa el marco de datos:

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

Finalmente, antes de retroceder la serie con esas características, se divide la serie en una participación de entrenamiento y prueba. Se establecen 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), ]

Modelling the series trend and seasonal components (Modelización de la tendencia de la serie y los componentes estacionales)

Primero se modela la tendencia de la serie haciendo una regresión de la serie con la variable de la 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 observa en el resultado de la 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, debido a que la mayor parte de la variación de la serie está relacionada con el patrón estacional.

Para situar en contexto, se usa el modelo creado para predecir los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La función de predicción predice los valores de una base de datos de entrada basada en un modelo dado.

Se crea una función de utilidad que represente la serie y la salida del modelo, utilizando el paquete plotly:

train$yhat <- predict(md_trend, newdata = train)

test$yhat <- predict(md_trend, newdata = test)

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

Se establecen 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 pudo 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, se mide la tasa de error del modelo tanto en los conjuntos de entrenamiento como 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 modelar y pronosticar el componente estacional sigue el mismo proceso que se aplica con la tendencia, realizando una regresión de la serie con la variable estacional que fue creada 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 que se realiza la regresión de la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías, las cuales son integradas con los valores de pendiente. Además, se nota que el R-cuadrado ajustado del modelo estacional es algo mayor con respecto al modelo de tendencia (0,78 en lugar de 0,1)

Antes de trazar el modelo y pronosticar los valores con la función plot_lm, se actualiza 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")

Como se evidencia en el gráfico anterior, el modelo realiza un trabajo positivo al capturar la estructura del patrón estacional de la serie. Sin embargo, se observa que falta la tendencia de la serie. Antes de agregar la tendencia como los componentes estacionales, se califica el desempeño 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 prueba está relacionada con el componente de tendencia que no fue incluido en el modelo. Seguidamente, se unen los dos componentes en un modelo y se pronostica los valores característicos de la serie:

md1 <- lm(y ~ seasonal + trend, data = train)

#Se mira 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, ocasiona un impulso adicional al R-cuadrado ajustado del modelo de 0,78 a 0,91. Lo anterior se observa 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 la puntuación MAPE del modelo tanto en las particiones de entrenamiento como 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 estacional generan un aumento significativo en la calidad del ajuste del modelo y en la precisión del modelo. No obstante, al observar la gráfica del ajuste y pronostico del modelo, se evidencia como la tendencia del modelo es demasiado lineal y falta 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 una mayor precisión del modelo.

Una técnica simple para capturar una tendencia no lineal es agregar un componente polinomial a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Se utiliza el argumento I, el cual permite aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Por ende, se usa 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 mejora significativa de la bondad de ajuste del modelo. En el otro modelo, como se observa 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 través 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")

Como se evidencia en el modelo que sigue la puntuación MAPE, la precisión del modelo obtuvo una mejora significativa al agregar la tendencia polinomial al modelo de regresión, debido a que el error en el conjunto de prueba 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

Hasta el momento se ha visto el proceso manual de convertir objetos ts en regresión lineal, Formato del modelo predictivo. La función tslm del paquete de pronóstico proporciona una función integrada que convierte el objeto ts en un modelo de predicción de regresión lineal. usar función tslm, puede configurar los componentes de regresión, así como otras funciones.

Ahora repetiremos el ejemplo anterior y predeciremos las últimas 12 observaciones de USgas 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$test

A continuación, aplicaremos la misma fórmula que usamos para crear el pronóstico md2 anterior

modelo usando 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 el resultado anterior, los dos modelos (md2 y md3) son idénticos.

Existen varias ventajas al usar la función tslm, en lugar de configurar 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 y ingeniería de características

-El objeto de salida soporta toda la funcionalidad del pronóstico (como elfunciones de precision y checkresiduals) y paquetes TSstudio (como elfunciones test_forecast y plot_forecast)

Modeling single events and non-seasonal events (Modelado de eventos únicos 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.

Dado que no se expresan en el modelo de regresión, tales eventos sesgarán los coeficientes de predicción ya que el modelo ponderará dichos eventos junto con los eventos regulares de la serie. El uso de variables codificadas, binarias o indicadoras puede ayudar al modelo a ignorar dichos eventos o ajustar los coeficientes del modelo en consecuencia.

Por ejemplo, en el gráfico de descomposición de la serie USgas que se muestra arriba, puede ver que la tendencia de la serie experimentó una ruptura estructural alrededor de 2010. Si bien el crecimiento fue relativamente plano antes de 2010, la pendiente de la tendencia cambió más tarde, con un crecimiento positivo. En este caso, podemos utilizar una variable binaria igual a cero para las observaciones antes de 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. El siguiente ejemplo ilustra el proceso de creación. de una variable binaria externa que vale 0 antes del año 2010 y 1 después, usando el Tabla USgas_df:

Como se puede observar en el resumen del modelo anterior, la variable quiebre estructural es estadísticamente significativa, con el nivel de 0,01. De manera similar, en el caso de un valor excepcional o un día festivo, se puede aplicar la codificación en caliente configurando una variable binaria en 1 siempre que ocurra un valor de excepción o ocurre un evento no estacional recurrente y 0 en caso contrario. Tenga en cuenta que una vez que haya entrenado el modelo de pronóstico con la función tslm usando variables externas, deberá generar los valores futuros de esas variables, ya que se usarán como entrada para el pronóstico.

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 evidenciar, 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.

Forecasting a series with multiseasonality components – a case study (Pronóstico de una serie con componentes de estacionalidad múltiple: 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 ofrece una amplia gama de opciones de personalización y nos permite modelar y pronosticar datos de series temporales complejas como las series multiestacionales.

En los siguientes ejemplos, utilizaremos la serie UKgrid para demostrar el enfoque de pronóstico de series multiestacionales 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 paquetes UKgrid. Esta serie representa una serie temporal de datos de alta frecuencia con una frecuencia de media hora. Usaremos la función extract_grid del paquete UKgrid para determinar la secuencia, características clave. (por ejemplo, formato de datos, variables, frecuencia, etc.) Esta función de conversión nos permite sumar la frecuencia de la serie de media hora a una frecuencia menor como horaria, diaria o mensual. Dado que nuestro objetivo aquí es pronosticar la demanda diaria para los próximos 365 días, estableceremos la serie en frecuencia diaria utilizando la estructura data.frame usamos:

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

Se observa patrones 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")

Preprocessing and feature engineering of the UK daily series (Preprocesamiento e ingeniería de funciones de la serie diaria del Reino Unido)

Para capturar los componentes estacionales de la serie, estableceremos la serie como frecuencia diaria y crearemos las siguientes dos características:

-Indicador del día de la semana

-Indicador del mes del año

Además, como es razonable suponer (confirmaremos esta suposición con la ACF función una vez que hemos 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: En los ejemplos a continuación, utilizaremos la serie UKgrid para ilustrar el enfoque de pronóstico de series multiestacionales con el modelo de regresión lineal.

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

Debido a que 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, además se usrá 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)

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 crear las nuevas funciones para el modelo y configurar el objeto ts, estamos listos para dividir la cadena de entrada que creamos (UKdaily) y el objeto de función externo correspondiente en una parte de entrenamiento y prueba. Dado que nuestro objetivo es predecir las próximas 365 observaciones, y la longitud de la serie es tan grande como (2540 observacione), podemos establecer la porción de prueba en la longitud del horizonte de predicción: 365 observaciones. Estableceremos la variable indicadora h en 365 y la usaremos para definir los segmentos y luego 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

De manera similar, necesitamos dividir las funciones que creamos para el modelo de regresión en (las características estacionales y de retraso) las particiones de entrenamiento y prueba exactamente en el mismo orden que usamos para el objeto ts correspondiente. Usaremos la función de índice data.frame para configurar la tabla UKdaily para entrenar y probar particiones:

train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]

Training and testing the forecasting model (Entrenamiento y prueba del modelo de pronóstico)

Una vez que hayamos creado las diversas funciones para el modelo, estamos listos para comenzar el proceso de capacitación 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)

En el gráfico de rendimiento anterior, podemos ver que el modelo de referencia está haciendo un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad de los días del año. De lo contrario, no puede captar las oscilaciones relativas a los días de la semana. Las puntuaciones MAPE del modelo, como podemos ver en el resultado de la siguiente función de precisión, son 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

Se intentara mejorar 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)

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 el gráfico de rendimiento del segundo modelo anterior, podemos ver la contribución de las características estacionales en el pronóstico, ya que el segundo modelo pudo capturar tanto las tendencias como los patrones multiestacionales en la serie. Esto también se puede ver 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

Se agrera la variable lag al modelo y se repite el mismo proceso:

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)

Se revisará el puntaje MAPE del tercer modelo, puesto que en la gráfica de rendimiento se ve diferencia significativa con respecto al segundo 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.

Model selection (Selección de modelo)

Ahora es el momento de elegir un modelo. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Dado que el segundo y el tercer modelo lograron puntajes MAPE bastante similares, con una pequeña ventaja para el tercer modelo, uno debe preguntarse si el uso de la demora variable vale la mejora de menos del 0,2 % en la tasa de error del conjunto de prueba. (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. También depende del objetivo del pronóstico. Se recomienda considerar la siguiente prueba:

La primera pregunta que se debe hacer en este caso es: ¿La variable rezagada es 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, 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 arrogado que la variable estadísticamente significativa

El tercer modelo puede ser más preciso por casualidad, ya que la diferencia es relativamente pequeña, no porque la variable de retraso contribuya a la precisión del modelo. Por lo tanto, la prueba retrospectiva de ambos modelos puede ayudar a confirmar si la contribución de la variable de retraso es consistente durante varios períodos de prueba. Os dejo que probéis ambos modelos como un divertido ejercicio.

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)

Residual analysis ( Análisis residual)

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 puede ver en el diagrama de resumen de residuos anterior, los residuos no son ruido blanco porque existe cierta autocorrelación entre la serie residual y sus retrasos. Técnicamente, esto es una indicación de que el modelo aún no ha capturado todos los patrones o la información que existe en la secuencia. Una forma de abordar este problema es identificar otras variables que puedan explicar el cambio en los residuos. El principal problema de este enfoque es que es difícil identificar variables externas que puedan explicar el cambio en los residuos y que además sean predecibles. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, pero es difícil predecir estos patrones climáticos hasta con un año de anticipación.

Un enfoque alternativo, cuando los patrones permanecen en los residuos del modelo, es tratar los residuos del modelo como datos de series de tiempo separados y modelarlos con un modelo de serie de tiempo predictivo diferente.

Finalizing the forecast (Finalizando el pronóstico)

Completemos el proceso y usemos el modelo entrenado seleccionado para predecir 365 observaciones futuras. Dado que estamos usando variables externas con la función tslm, necesitaremos restar sus valores futuros. Es relativamente simple ya que usamos 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 simple es generar primero las fechas correspondientes de las observaciones del pronóstico y luego restar el día de la semana y el mes del año de este objeto:

UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1), 
                                       by = "day", 
                                       length.out = h))

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 este capítulo, presentamos las aplicaciones predictivas del modelo de regresión lineal. Aunque el modelo de regresión lineal no está diseñado para manejar datos de series de tiempo, la ingeniería de características simple nos permite convertir un problema de pronóstico en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal sobre otros modelos tradicionales de series temporales es la capacidad de incorporar variables y factores extrínsecos del modelo.

En este capítulo, presentaremos aplicaciones de pronóstico del modelo de regresión lineal. Aunque el modelo de regresión lineal no fue diseñado para trabajar con datos de series de tiempo, la ingeniería de características simple nos permite transformar el problema de pronóstico en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal sobre otros modelos tradicionales de series temporales es la capacidad del modelo para incorporar variables y factores externos.