Sección 9
Angello Restrepo
2023-04-17
0.1 Previsiones con regresión lineal
El modelo de regresión lineal es uno de los métodos más comunes para identificar y cuantificar la relación entre una variable dependiente y una única (regresión lineal univariante) o múltiples (regresión lineal multivariante) variables independientes. Este modelo tiene una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo y, en particular, la previsión de series temporales. En este capítulo se tratan los siguientes temas - Enfoques de previsión con modelos de regresión lineal - Extracción y estimación de los componentes de la serie - Tratamiento de rupturas estructurales, valores atípicos y eventos especiales - Previsión de series con multiestacionalidad ###Regresión lineal El modelo expresa la variable dependiente como una combinación lineal de las variables independientes. Una relación lineal entre las variables dependiente e independiente puede generalizarse mediante las siguientes ecuaciones.
En el caso de una única variable independiente, la ecuación es la siguiente: \(Y_i = \beta_0 + \beta_1 * X_i + \epsilon_i\). En el caso de un numero n de variables independientes, la ecuacion es: \(Y_i = \beta_0 + \beta_1 X_1,_i + \beta_2 X_2,_i +...+ \beta_n X_n,_i +\epsilon_i\).
Suponiendo que las ecuaciones anteriores representen la verdadera naturaleza de la relación lineal entre las variables dependiente e independiente, el modelo de regresión lineal proporciona una estimación de dichos coeficientes (es decir, ), que puede formalizarse mediante las siguientes ecuaciones:
- Para el modelo de regresión lineal univariante, la ecuación es la siguiente: \(\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 X_i + \epsilon_i\). -Para el modelo de regresión lineal multivariante, la ecuación es la siguiente: \(\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 X_i + \hat{\beta}_2 X_2,_i +...+ \hat{\beta}_n X_n,_i + \epsilon_i\).
La estimación de los coeficientes del modelo se basa en los dos pasos siguientes: - Definir una función de coste (también conocida como función de pérdida): establecer una métrica de error para minimizarla. - Aplicar la optimización matemática para minimizar la función de coste. El método de estimación más común consiste en aplicar el método de mínimos cuadrados ordinarios (OLS) como técnica de optimización para minimizar la suma de los cuadrados de los residuos ( ).
como técnica de optimización para minimizar la suma de cuadrados de los residuos ( ). Elevar los residuos al cuadrado tiene dos efectos: - Evita que los valores positivos y negativos se anulen al sumarlos. - El efecto cuadrado proporciona una penalización exponencial para los residuos con mayor distancia, ya que su coste es mayor. ###Estimación de coeficientes con el método de minimos cuadrados ordinarios(OLS) El método de mínimos cuadrados ordinarios (MCO) se utiliza para identificar los coeficientes que minimizan la suma de cuadrados residuales de un modelo de regresión lineal. El método MCO utiliza álgebra lineal básica y cálculo para resolver los coeficientes. Para ello, la función de coste se define como la suma de las diferencias al cuadrado entre los valores reales y los predichos. El método MCO se aplica para transformar la función de coste en forma de matriz para simplificar. Al diferenciar la función de coste con respecto a los coeficientes y ponerla a cero, se pueden resolver los coeficientes que minimizan la función de coste. Las propiedades clave de MCO son la insesgadez de la estimación de coeficientes, que la línea de regresión de la muestra pase por la media de X e Y y que la media de los residuos sea igual a cero.
0.1.1 Los supuestos del método de minimos cuadrados ordinarios(OLS)
El modelo de mínimos cuadrados ordinarios (MCO) tiene ciertos supuestos que deben cumplirse para que el modelo sea válido. Estos supuestos incluyen la linealidad de los coeficientes del modelo, la ausencia de colinealidad perfecta entre las variables independientes, la varianza distinta de cero de todas las variables independientes y que el término de error sea independiente e idénticamente distribuido (i.i.d) con media O y varianza constante. Además, las variables dependientes e independientes deben extraerse de la población en una muestra aleatoria. Estos supuestos pueden incumplirse en algunos casos, como con los datos de series temporales en los que puede existir cierto grado de correlación en los residuos del modelo. En tales casos, pueden utilizarse métodos alternativos, como los modelos ARIMA.
0.1.2 Previsiones con regresión lineal
El modelo de regresión lineal es un modelo genérico utilizado para una amplia gama de aplicaciones, incluidos el análisis predictivo y la inferencia causal. No se diseñó explícitamente para datos de series temporales, a diferencia de los modelos tradicionales de series temporales como ARIMA o Holt-Winters. Para utilizar la regresión lineal en la previsión de series temporales, el primer paso consiste en identificar las características clave, los patrones, los valores atípicos y otras características de la serie. El segundo paso consiste en transformar estas características en variables de entrada y utilizar la regresión para crear un modelo de previsión. Las principales características de un modelo de previsión de regresión lineal son los componentes de tendencia y estacionales, que deben identificarse y transformarse en variables de entrada.
0.1.3 Previsión de la tendencia y los componentes estacionales
En el análisis de datos de series temporales, es esencial entender los componentes estructurales que conforman la serie, a saber, los componentes de tendencia, estacionales, de ciclo e irregulares. En el capítulo anterior, se presentó la función descomponer y se demostró cómo se pueden descomponer los componentes estructurales de una serie utilizando esta función. Además, se proporcionaron ecuaciones para la descomposición de una serie con una estructura aditiva o multiplicativa.
Sin embargo, el análisis de la estructura de una serie no es suficiente para hacer predicciones precisas. Es aquí donde entra en juego la regresión lineal, que permite crear un modelo de predicción utilizando los componentes estructurales como variables independientes y el término de error como el componente irregular. En este capítulo, se describe cómo se pueden transformar las ecuaciones de descomposición en un modelo de regresión lineal, con los componentes de tendencia y estacionales como variables independientes y el término de error que representa el componente irregular.
Además, se discute la transformación de una serie con una estructura multiplicativa en una estructura aditiva utilizando la transformación logarítmica. En general, el pasaje se centra en la importancia de entender los componentes estructurales de una serie temporal y cómo se pueden utilizar para crear un modelo de predicción preciso.
0.1.4 Características técnicas de los componentes de serie
Para crear las entradas de regresión, se debe entender la estructura de la serie, lo cual se ejemplificará con la serie del gas estadounidense. Carguemos de nuevo la serie desde el paquete TSstudio y representémosla con ts_plot :
library(TSstudio)
data("USgas")
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Además, revisemos 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 10La serie del gas de EE.UU. es mensual y tiene un componente estacional fuerte y una tendencia estable. Podemos analizar su estructura con ts_decompose.
ts_decompose(USgas)La tendencia de la serie de tiempo no es estrictamente lineal entre 2000 y 2018, lo que es importante al crear la entrada de tendencia para el modelo de regresión.
Antes de utilizar la función lm, la función de regresión lineal incorporada en R del paquete stats, tendremos que transformar la serie de un objeto ts a un objeto data. frame. Por lo tanto, utilizaremos la función ts_to_prophet del paquete TSstudio:
USgas_df <- ts_to_prophet(USgas)La función transforma el objeto ts en dos columnas de data. frame, donde las dos columnas representan los componentes temporal y numérico de la serie, respectivamente:
head(USgas_df)Se crea la entrada de la regresión para el modelo transformando la serie en un data frame. La primera característica que se crea es la tendencia de la serie mediante la indexación cronológica de las observaciones.
USgas_df$trend <- 1:nrow(USgas_df)Usando el índice de la serie, podemos estimar el crecimiento marginal de un mes a otro. Luego, para crear la característica de componente estacional, se utiliza una variable categórica con 12 categorías, correspondientes a cada mes del año, que se extraen de la variable de fecha usando la función month del paquete lubridate.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionUSgas_df$seasonal <- month(USgas_df$ds, label = T)Revisemos ahora el marco de datos después de añadir las nuevas funciones:
head(USgas_df)Antes de empezar a realizar la regresión de la serie con esas características, dividiremos la serie en una partición de entrenamiento y otra de prueba. Estableceremos los últimos 12 meses de la serie como 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), ]Revisaremos cómo el modelo de regresión captura los componentes por separado y juntos, tras crear los marcos de datos de entrenamiento y prueba.
0.1.5 Modelización de la tendencia de la serie y los componentes estacionales
Primero modelizaremos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, sobre la partición de entrenamiento y utilizaremos la función summary para revisar los detalles del modelo:
md_trend <- lm(y ~ trend, data = train)
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-09Se usará la función “summary” para revisar los detalles del modelo de regresión. El coeficiente de la variable de tendencia es significativo, pero el R-cuadrado ajustado es bajo debido al patrón estacional dominante en la serie. El valor p indica la probabilidad de rechazar la hipótesis nula. Se recomienda visualizar los datos a través de la predicción de valores ajustados y pronosticados con la función “predict”.
Lo utilizaremos para predecir los valores ajustados y previstos del modelo de tendencia que hemos entrenado antes:
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)Vamos a crear una función de utilidad que traza la serie y la salida del modelo, utilizando el paquete plotly:
library(ggplot2)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = ""),
yaxis = list(title = "Billion Cubic Feet"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Los argumentos de la función son los siguientes: - datos: Los datos de entrada, un objeto data. frame siguiendo la misma estructura que la del USgas_df (incluyendo la variable yhat) - train: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo - test: Del mismo modo, el conjunto de prueba correspondiente que se utilizó para evaluar el modelo de previsión - title: El título del gráfico, por defecto NULL.
Establezcamos 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 capturó la tendencia general pero no la ruptura estructural de la tendencia en 2010. Se presentará un método avanzado para capturar una tendencia no lineal.
Por último, para el análisis comparativo, queremos medir la tasa de error del modelo tanto en el conjunto de entrenamiento como en el de pruebas:
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.1299951Modelar el componente estacional implica regresar la serie utilizando la variable categórica creada anteriormente, tal como se hizo con la tendencia. Repasemos los detalles del modelo:
md_seasonal <- lm(y ~ seasonal, data = train)
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-16Dado que regresionamos 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 que están incrustadas con los valores de la pendiente.
Antes de trazar el modelo ajustado y los valores previstos con la función plot_lm, actualizaremos los valores de yhat con la función predict. Luego podemos utilizar la función plot_lm para visualizar el modelo ajustado y los valores previstos:
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 puede verse en el gráfico anterior, el modelo capta bastante bien la estructura del patrón estacional de la serie. Sin embargo, se puede observar que falta la tendencia de la serie. Antes de añadir tanto la tendencia como los componentes estacionales, puntuaremos 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.21502100El error en la prueba se debe a la falta de la componente de tendencia en el modelo. Ahora se procederá a combinar ambos componentes en un modelo y hacer pronósticos. Luego repasemos el resumen del modelo tras realizar una regresión de la serie con los componentes de tendencia y estacional:
md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -514.73 -77.17 -17.70 85.80 336.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1733.7153 17.0794 101.509 < 2e-16 ***
## seasonal.L -498.1709 29.6354 -16.810 < 2e-16 ***
## seasonal.Q 1115.2951 29.5872 37.695 < 2e-16 ***
## seasonal.C 78.9835 29.5103 2.676 0.00802 **
## seasonal^4 175.6505 29.5196 5.950 1.09e-08 ***
## seasonal^5 285.0192 29.5568 9.643 < 2e-16 ***
## seasonal^6 -49.3611 29.5319 -1.671 0.09610 .
## seasonal^7 -192.3050 29.4540 -6.529 4.77e-10 ***
## seasonal^8 81.8787 29.3835 2.787 0.00581 **
## seasonal^9 44.4849 29.3480 1.516 0.13106
## seasonal^10 76.8636 29.3372 2.620 0.00943 **
## seasonal^11 -11.2755 29.3353 -0.384 0.70109
## trend 2.6182 0.1305 20.065 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared: 0.9142, Adjusted R-squared: 0.9094
## F-statistic: 189.2 on 12 and 213 DF, p-value: < 2.2e-16La regresión de la serie con la tendencia y los componentes estacionales juntos proporciona un aumento adicional de la R-cuadrado ajustada del modelo de O . 78 a O . 91.Esto puede observarse en el gráfico de los resultados 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")Midamos la puntuación MAPE del modelo tanto en la partición de entrenamiento como en la 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.09143290Mejora la regresión de la serie al agregar la tendencia y los componentes estacionales, pero aún falta capturar la ruptura estructural de la tendencia. Agregar un componente polinómico puede mejorar la precisión. Se puede utilizar el argumento I para añadir un segundo grado del polinomio a la tendencia y capturar su curvatura a lo largo del tiempo.
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)
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-16El polinomio de segundo grado no mejoró significativamente el modelo. Sin embargo, en el otro modelo, se pudo capturar la ruptura estructural de la tendencia con este cambio.
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 podemos ver en el modelo que sigue a la puntuación MAPE, la precisión del modelo mejoró significativamente al añadir la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas bajó 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.042126180.1.6 La función tslm
Utilizando la función tslm del paquete forecast, podemos transformar un objeto ts en un modelo de previsión de regresión lineal. Esto nos permite establecer el componente de regresión y otras características. Utilizaremos la función tslm para pronosticar las 12 últimas observaciones de la serie USgas, utilizando los componentes de tendencia, tendencia al cuadrado y estacional. Primero dividimos la serie en conjuntos de entrenamiento y de prueba utilizando la función ts_split.
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testA continuación, aplicaremos la misma fórmula que utilizamos para crear el modelo de previsión md2 anterior utilizando la función ts lm:
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoomd3 <- tslm(train.ts ~ season + trend + I(trend^2))
summary(md3)##
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## season2 -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## season3 -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## season4 -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## season5 -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## season6 -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## season7 -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## season8 -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## season9 -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## season10 -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## season11 -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## season12 -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16La función tslm es una forma eficiente y fácil de transformar un objeto ts en un modelo de previsión de regresión lineal, ya que no requiere la transformación de la serie en un objeto data.frame. Además, el objeto de salida admite toda la funcionalidad de los paquetes forecast y TSstudio. En comparación con establecer manualmente un modelo de regresión, la función tslm ofrece ventajas como eficiencia y menor necesidad de ingeniería de características.
0.1.7 Modelización de acontecimientos aislados y no estacionales
En los datos de series temporales, pueden haber patrones inusuales como valores atípicos, rupturas estructurales y eventos recurrentes no estacionales que sesgarán los coeficientes del modelo si no se expresan en él. Se puede utilizar la codificación en caliente con variables binarias para ajustar los coeficientes del modelo o ignorar estos eventos. Por ejemplo, una variable binaria se puede utilizar para capturar la ruptura estructural de la tendencia en la serie USgas alrededor del año 2010.
Para regresar un modelo tslm con variables externas, se necesita un objeto data.frame con las variables correspondientes. Se puede crear una variable binaria externa que sea igual a 0 antes del año 2010 y 1 después. El ejemplo muestra cómo hacerlo con la tabla USgas_df y luego remodelar la serie USgas. La función de resumen se puede utilizar para revisar la salida del modelo.
r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse(year(USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)
summary(md3)##
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break,
## data = USgas_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -469.25 -50.68 -2.66 63.63 275.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.661e+03 3.200e+01 83.164 < 2e-16 ***
## season2 -3.054e+02 3.448e+01 -8.858 2.61e-16 ***
## season3 -4.849e+02 3.448e+01 -14.062 < 2e-16 ***
## season4 -9.272e+02 3.449e+01 -26.885 < 2e-16 ***
## season5 -1.108e+03 3.449e+01 -32.114 < 2e-16 ***
## season6 -1.127e+03 3.450e+01 -32.660 < 2e-16 ***
## season7 -9.568e+02 3.450e+01 -27.730 < 2e-16 ***
## season8 -9.340e+02 3.451e+01 -27.061 < 2e-16 ***
## season9 -1.138e+03 3.452e+01 -32.972 < 2e-16 ***
## season10 -1.040e+03 3.453e+01 -30.122 < 2e-16 ***
## season11 -7.896e+02 3.497e+01 -22.577 < 2e-16 ***
## season12 -2.649e+02 3.498e+01 -7.571 9.72e-13 ***
## trend -1.928e+00 4.479e-01 -4.304 2.51e-05 ***
## I(trend^2) 1.862e-02 1.676e-03 11.113 < 2e-16 ***
## s_break 6.060e+01 2.836e+01 2.137 0.0337 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109 on 223 degrees of freedom
## Multiple R-squared: 0.9423, Adjusted R-squared: 0.9387
## F-statistic: 260.3 on 14 and 223 DF, p-value: < 2.2e-16Resumen: En el modelo de regresión lineal tslm, se pueden utilizar variables externas para tener en cuenta patrones inusuales en los datos, como rupturas estructurales, valores atípicos o días festivos. Esto se puede lograr mediante la codificación en caliente de variables binarias. Después de entrenar el modelo, es necesario tener en cuenta las variables externas al producir predicciones futuras.
0.1.8 Previsión 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 ofrece una amplia gama de opciones de personalización y nos permite modelizar y pronosticar datos de series temporales complejas, como las series con multiestacionalidad. En los siguientes ejemplos, utilizaremos la serie UKgrid para demostrar el enfoque de previsión de una serie multiestacional con un modelo de regresión lineal.
0.1.9 Las series UKgrid
La serie UKgrid representa la demanda de electricidad en el Reino Unido y se puede transformar con la función extract_grid del paquete UKgrid para definir sus características y frecuencia. Se establecerá la serie en frecuencia diaria para prever la demanda diaria en los próximos 365 días. Usaremos la función head para revisar las variables de la serie.
library(UKgrid)
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
head(UKdaily)Como puedes ver, esta serie tiene dos variables: - TIMES TAMP: Un objeto de fecha utilizado como marca de tiempo o índice de la serie - ND: La demanda neta de electricidad
Utilizaremos la función ts_plot para trazar y revisar la estructura de la serie:
ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")Como se puede ver en el gráfico anterior, la serie tiene una clara tendencia a la baja y tiene un patrón estacional de cadena. Esta serie tiene múltiples patrones de estacionalidad: - Diaria: 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 puede verse en el siguiente heatrnap 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")Como puede ver en el mapa térmico de la serie, la demanda global aumenta a lo largo de las semanas de invierno (por ejemplo, las semanas naturales 1 a 12 y 44 a 52). Además, se puede observar el cambio de la serie durante los días de las semanas, ya que la demanda aumenta durante los días laborables de la semana y disminuye durante el fin de semana.
0.1.10 Ingenieria de preprocesamiento y caracteristicas de la serie UKdaily
Crear características como un indicador del día de la semana y del mes del año para capturar los componentes estacionales de la serie, y también una variable de rezago de 365 observaciones para capturar la fuerte correlación con los rezagos estacionales. Se utilizará el paquete dplyr para crear estas características.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionUKdaily <- UKdaily %>%
mutate(wday = wday(TIMESTAMP, label = TRUE),
month = month(TIMESTAMP, label = TRUE),
lag365 = dplyr::lag(ND, 365)) %>%
filter(!is.na(lag365)) %>%
arrange(TIMESTAMP)
str(UKdaily)## 'data.frame': 4939 obs. of 5 variables:
## $ TIMESTAMP: Date, format: "2006-04-01" "2006-04-02" ...
## $ ND : int 1718405 1691341 1960325 2023886 2026204 2008422 1981175 1770440 1749715 2012865 ...
## $ wday : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 7 1 2 3 4 5 6 7 1 2 ...
## $ month : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 4 4 4 4 4 4 4 4 4 4 ...
## $ lag365 : int 1920069 1674699 1631352 1916693 1952082 1964584 1990895 2003982 1811436 1684720 ...Como la entrada de la función tslm debe estar en formato ts, convertiremos la serie en un objeto ts. Utilizaremos la primera marca de tiempo de la serie y las funciones year y day (el día del año) del paquete lubridate para establecer el punto de inicio del objeto:
start_date <- min(UKdaily$TIMESTAMP)
start <- c(year(start_date), yday(start_date))Después de transformar la serie en un objeto ts, podemos volver atrás y confirmar la suposición que hicimos sobre el nivel de correlación entre la serie y sus rezagos estacionales con la función ts_acf (una versión interactiva de la función acf del paquete TSstudio). Revisaremos la correlación de la serie con sus retardos de los últimos cuatro años:
UK_ts <- ts(UKdaily$ND,
start = start,
frequency = 365)
ts_cor(UK_ts, lag.max = 365 * 4)Dividimos la serie de entrada y el objeto de características externas en una partición de entrenamiento y otra de prueba, estableciendo la partición de prueba en la longitud del horizonte de previsión (365 observaciones). La función ts_split se utiliza para dividir la serie y las características.
h <- 365
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$test
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]###Entrenamiento y prueba del modelo de previsión Una vez creadas las distintas características del modelo, estamos listos para iniciar el proceso de entrenamiento del modelo de previsión. - Modelo base: Regresión de la serie con el componente estacional y de tendencia utilizando las características incorporadas de la función tslm. - Modelo multiestacional con desfase estacional: Utilizando, además de los indicadores estacionales, la variable de desfase estacional. - Rendimiento del modelo en el conjunto de entrenamiento y de pruebas mediante la puntuación MAPE. Empezar con un modelo simplista (modelo de referencia) nos permitirá observar si añadir nuevas características contribuye al rendimiento del modelo o si debemos evitarlo, ya que añadir más características o complejidad al modelo no siempre produce mejores resultados.
md_tslm1 <- tslm(train_ts ~ season + trend)
fc_tslm1 <- forecast(md_tslm1, h = h)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)El gráfico muestra que el modelo de referencia captura bien la tendencia y estacionalidad diaria, pero no la oscilación semanal. El MAPE del modelo es 6,09% en la partición de entrenamiento y 7,77% en la de prueba, según la función de precisión.
accuracy(fc_tslm1, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -4.781286e-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.071899Intentaremos mejorar el modelo agregando las características del día de la semana y el mes del año. Para hacer esto, especificamos los datos de entrada con el argumento ‘datos’ y repetimos el proceso anterior de entrenar el modelo y evaluar su rendimiento con una previsión.
md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)
fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
forecast.obj= fc_tslm2, test = test_ts)El segundo modelo mejora la precisión al incluir características estacionales y se refleja en el gráfico y la puntuación MAPE (2,87% y 5,23% en entrenamiento y prueba, respectivamente).
accuracy(fc_tslm2, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set 8.172823e-12 70245.98 52146.79 -0.1738708 3.167605 0.4370853
## Test set -1.764563e+04 80711.71 65373.21 -1.3715505 4.682071 0.5479470
## ACF1 Theil's U
## Training set 0.7513664 NA
## Test set 0.6075598 0.68445Por último, pero no menos importante, vamos a añadir la variable lag al modelo, y repetir el mismo proceso que antes. El rendimiento del tercer modelo puede verse en el siguiente gráfico:
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 en el gráfico de rendimiento del tercer modelo si existe una diferencia significativa con respecto al segundo. Por lo tanto, revisemos la puntuación MAPE del tercer modelo:
accuracy(fc_tslm3, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.836939e-13 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
## ACF1 Theil's U
## Training set 0.7500146 NA
## Test set 0.6094666 0.6925767Los resultados del tercer modelo muestran una pequeña mejora en la precisión del modelo, con un 2,81% en el conjunto de entrenamiento y un 5,07% en el conjunto de pruebas.
0.1.11 Selección del modelo
Se deben elegir entre el segundo y el tercer modelo, ya que son superiores al primero en términos de precisión. La elección entre los dos modelos dependerá del costo adicional y de la importancia de una pequeña mejora en la tasa de error. Se recomienda evaluar si la variable de retardo es estadísticamente significativa mediante la función de resumen.
summary(md_tslm3)$coefficients %>% tail(1)## Estimate Std. Error t value Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05El valor p de la variable lag365 indica 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 y comprobar si la variable adicional es significativa
anova(md_tslm3)Seleccionamos el tercer modelo debido a su mayor precisión, pero podría aplicarse un enfoque más robusto de selección de características. Ahora, debemos entrenar el modelo en toda la serie y prever los próximos 365 días.
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365,
data = UKdaily)0.1.12 Análisis de residuos
Realizamos 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-16Los residuos no son ruido blanco debido a la autocorrelación, lo que indica que el modelo no capturó toda la información en las series. Identificar variables adicionales que puedan explicar la variación de los residuos es difícil. Una alternativa es tratar los residuos como datos independientes y modelizarlos con otro modelo de previsión de series temporales, como una regresión con error ARIMA en el Capítulo 11.
0.1.13 Finalizando la previsión
Para pronosticar las 365 observaciones futuras utilizando el modelo seleccionado, necesitamos generar los valores futuros de las variables externas utilizadas en el modelo. Esto es fácil de hacer ya que utilizamos variables deterministas. Creamos un data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras y extraemos el día de la semana y el mes del año de las fechas correspondientes a las observaciones pronosticadas.
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1),
by = "day",
length.out = h))A continuación, podemos utilizar la variable date para crear las variables wday y month con el paquete paquete lubridate:
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)Utilicemos la función forecast para crear la previsión:
UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)Por último, trazaremos la previsión final con la función plot_forecast del paquete TSstudio:
plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")0.1.14 Resumen
En este capítulo hemos presentado las aplicaciones de previsión del modelo de regresión lineal. Aunque el modelo de regresión lineal no se diseñó para manejar datos de series temporales, con una simple ingeniería de características podemos transformar un problema de previsión en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal con respecto a otros modelos tradicionales de series temporales es su capacidad para incorporar variables y factores externos. No obstante, este modelo puede manejar series temporales con patrones multiestacionales, como vimos con la previsión de la demanda de electricidad en el Reino Unido. Por último, pero no por ello menos importante, los enfoques de previsión que hemos demostrado en este capítulo serán la base para la modelización avanzada con modelos de aprendizaje automático que trataremos en el Capítulo 12, Previsión con modelos de aprendizaje automático. En el próximo capítulo, presentaremos los métodos de suavización exponencial, una familia de modelos de previsión basados en un enfoque de media móvil ponderada.