El modelo de regresión linear es uno de los más comunes para cuantificar e identificar las relaciones que existen entre una variable dependiente con una o varias variables independientes. Por lo que, este modelo expresa una varible dependiente y su relación linear (de sus coeficientes) respecto a una o varias independientes.
Para estimar los coeficientes de este modelo se debe definir una función (cost function) y luego minimizarla mediante la optimización matemática. Para realizar esto último, el método más común es el de los mínimos cuadrados lineales (Ordinary Least Squares), el cual está basado en el álgebra linear y cálcula, e identifica los coeficientes para así minimizar la suma residual de cuadrados. Esto tiene dos efectos, evitar que los valores positivos y negativos se cancelen entre sí al sumarlos y a su vez aumentar el valor de los residuos con mayor distancia.
Las propiedades clave del modelo de mínimos cuadrados lineales son las siguientes: -La insesgadez de la estimación de los coeficientes con respecto a los valores reales. -La línea de regresión de la muestra siempre pasará por la media de X y Y. -La media de la estimación hecha con los mínimos cuadrados lineales de la variable dependiente, es igual a la media de la variable dependiente. -La media del vector de residuos es equivalente a cero.
-Los coeficientes deben seguir una estructura linear. -Ninguna de las variables independientes es una combinación lineal de ninguna de las otras variables independientes. -Todas las variables independientes deben tener una varianza no constante. -Tanto la variable dependiente como la independiente se extraen de la población en una muestra aleatoria.
La regresión linear es un modelo genérico con una amplia gama de aplicaciones, por lo que, no fue diseñado explícitamente para manejar y pronosticar datos de series temporales.
La previsión identifica la estructura de la serie, los patrones, los valores atípicos y otras características. Al ver esas características en variables de entrada y hacer una regresión con las series podemos crear un modelo de previsión.
Así, las características principales de un modelo de previsión de regresión lineal son la tendencia y los componentes estacionales.
Es importante la identificación de los componentes de tendencia y estacionales de las series y su posterior transformación en variables de entrada del modelo de regresión.
La tendencia (Trend), representa el crecimiento de la serie respecto al tiempo tras ajustar y eliminar los efectos estacionales. Respecto a lo estacional (Seasonal), esto implica un patrón cíclico recurrente que se deriva de las unidades de frecuencia de la serie. El ciclo (Cycle) se refiere a un patrón cíclico que no está relacionado con la unidad de frecuencia de la serie. Lo irregular, se refiere a cualquiera otra pauta que no sea captada por los componentes de tendencia, estacional y componentes del ciclo.
Sin embargo, podemos dejar de lado el componente de ciclo, ya que normalmente se fusiona con el componente de tendencia (o se ignora el componente de ciclo).
La transformación de una serie con estructura multiplicativa en una formación de regresión lineal requiere una transformación de la serie en una estructura aditiva. Una vez transformada la serie en una estructura aditiva, la transformación a una formación de regresión es sencilla.
Después de convertir la serie a una estructura aditiva, la conversión a la formación de regresión lineal es sencilla y sigue el mismo proceso que el anterior.
library(TSstudio)
data("USgas")
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Luego de esto, se representan las principales características de la serie.
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 10ts_decompose(USgas)USgas_df <- ts_to_prophet(USgas)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.1Luego de este proceso se determina la tendencia (Trend) de la serie
USgas_df$trend <- 1:nrow(USgas_df)Posterior a esto, se hace la parte estacional.
head(USgas_df)## ds y trend
## 1 2000-01-01 2510.5 1
## 2 2000-02-01 2330.7 2
## 3 2000-03-01 2050.6 3
## 4 2000-04-01 1783.3 4
## 5 2000-05-01 1632.9 5
## 6 2000-06-01 1513.1 6h <- 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), ]En primer lugar, se modela la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia:
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-09El coeficiente de la variable de tendencia es estadísticamente significativo a un nivel de 0,001, ya que la mayor parte de la variación de la serie está relacionada con el patrón estacional.
La función de predicción del paquete stats, como su nombre indica, predice los valores de una base de datos de entrada basándose en un modelo determinado.
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)Modelado de salida:
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = ""),
yaxis = list(title = "Billion Cubic Feet"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Entradas de la función con la salida del modelo.
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")Se mide la tasa de error del modelo, 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.1299951La modelización y predicción del componente estacional sigue el mismo procedimiento que utilizamos para la tendencia.
La función tslm del paquete de previsión proporciona una función incorporada para transformar un objeto ts en un modelo de previsión de regresión lineal. Con el ejemplo anterior se prevee las últimas 12 observaciones de la serie estadounidense con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales.
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testlibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoomd3 <- tslm(train.ts ~ season + trend + I(trend^2))
summary(md3)##
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -468.47 -54.66 -2.21 63.11 294.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## season2 -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## season3 -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## season4 -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## season5 -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## season6 -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## season7 -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## season8 -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## season9 -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## season10 -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## season11 -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## season12 -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
## trend -1.270e+00 4.472e-01 -2.840 0.00494 **
## I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
## F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16Como puede verse en el resultado anterior, los dos modelos son idénticos. Por lo que, el uso de la función tslm tiene varias ventajas en lugar de crear manualmente un modelo de regresión para las series con la función lm: En primer lugar, no es necesario transformar las series en un objeto data.frame y crear funciones.Seguidamente, la salida es compatible con todas las funcionalidades de previsión (por ejemplo, las funciones accuracy y checkresiduals) y los paquetes de TSstudio (por ejemplo, las funciones test_forecast y plot_forecast).
Los datos de las series temporales pueden contener patrones inusuales que se repiten o no en el tiempo.
La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejercicio muestra el proceso de creación de una variable binaria externa igual a 0 antes de 2010 y 1 después de 2010 utilizando la tabla USgas_df.
r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse((USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1
#Se usa la siguiiente linea de código para remodelar la serie USgas
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)
#A continuación se usa la función summary para el modelo output
summary(md3)##
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break,
## data = USgas_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -473.54 -55.44 -0.79 63.31 285.69
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.648e+03 3.166e+01 83.653 < 2e-16 ***
## season2 -3.058e+02 3.476e+01 -8.799 3.77e-16 ***
## season3 -4.857e+02 3.476e+01 -13.974 < 2e-16 ***
## season4 -9.283e+02 3.476e+01 -26.709 < 2e-16 ***
## season5 -1.109e+03 3.476e+01 -31.912 < 2e-16 ***
## season6 -1.129e+03 3.476e+01 -32.469 < 2e-16 ***
## season7 -9.591e+02 3.476e+01 -27.591 < 2e-16 ***
## season8 -9.366e+02 3.476e+01 -26.943 < 2e-16 ***
## season9 -1.141e+03 3.477e+01 -32.829 < 2e-16 ***
## season10 -1.044e+03 3.477e+01 -30.017 < 2e-16 ***
## season11 -7.927e+02 3.522e+01 -22.507 < 2e-16 ***
## season12 -2.683e+02 3.522e+01 -7.618 7.18e-13 ***
## trend -1.560e+00 4.168e-01 -3.743 0.000231 ***
## I(trend^2) 1.869e-02 1.689e-03 11.063 < 2e-16 ***
## s_break NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 109.9 on 224 degrees of freedom
## Multiple R-squared: 0.9412, Adjusted R-squared: 0.9377
## F-statistic: 275.6 on 13 and 224 DF, p-value: < 2.2e-16Como se puede ver en el resumen del modelo anterior, la variable de rotura estructural es estadísticamente significativa, con un nivel de 0,01. Del mismo modo, en el caso de valores atípicos o días festivos, la codificación en caliente se puede aplicar estableciendo una variable binaria que sea igual a 1 cada vez que ocurra un evento atípico o no estacional, y 0 en caso contrario.
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 representa la demanda nacional de electricidad en el Reino Unido. 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. Puesto que nos permite agregar la frecuencia de la serie de media hora a una frecuencia más baja, como horaria, diaria o mensual.
library(UKgrid)
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
head(UKdaily)## TIMESTAMP ND
## 1 2005-04-01 1920069
## 2 2005-04-02 1674699
## 3 2005-04-03 1631352
## 4 2005-04-04 1916693
## 5 2005-04-05 1952082
## 6 2005-04-06 1964584ts_plot(UKdaily, title = "The UK National Demand for Electricity", Ytitle = "MW", Xtitle = "Year")Los patrones diarios, semanales y mensuales se pueden observar mediante la función ts_heatmap del paquete TSstudio. ## Pre-procesamiento e ingeniería de características de la serie UKdaily
Para poder observar los datos estacionales de la serie se necesita establecer la serie como frecuencia diaria, junto con dos indicadores; dia de la semana y mes del año.
Finalmente, se dividen las características que se crearon para el modelo de regresión en una partición de entrenamiento y prueba.
Después de haber creado las diferentes características para el modelo, se puede comenzar el proceso de entrenamiento del modelo de pronóstico. Utilizaremos la partición de entrenamiento y comenzaremos a entrenar los siguientes tres modelos:
Modelo de línea base: Regresión de la serie con el componente estacional y de tendencia utilizando las características integradas de la función tslm.
Modelo multiestacional: Se suman 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 los indicadores estacionales y la variable de rezago estacional.
Es difícil observar si hay una diferencia significativa del tercer modelo con respecto al segundo modelo. Para elegir el indicado, no hay un método correcto, ya que depende del objetivo del pronóstico. Se recomienda que considere la siguiente pregunta ¿Es la variable de retraso estadísticamente significativa?
Si esta no es estadísticamente significativa, no tiene sentido continuar, por lo que es mejor abandonar la variable. En el caso del tercer modelo, podemos utilizar la función de resumen para observar el nivel de significancia de la variable lag365. De la misma manera, podemos aplicar una prueba ANOVA con la función anova del paquete stats, y comprobar si la variable adicional es significativa.
Sin embargo, es posible que el tercer modelo sea más preciso sólo por casualidad, y no porque las variables delay contribuyan a la precisión del modelo, ya que las diferencias son relativamente pequeñas. Por lo tanto, el backtesting de los dos modelos puede ayudar a verificar que la contribución de las variables delay es coherente.
Justo antes de finalizar el pronóstico, se debe realizar un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals.
La existencia de cierta autocorrelación entre las series residuales y sus retrasos, es técnicamente una indicación de que el modelo no capturó todos los patrones o información que existe en la serie. Para abordar este problema se debe identificar variables adicionales que puedan explicar la variación en los residuos. Como enfoque alternativo, cuando los patrones sobrantes en los residuos del modelo se debe tratar los residuos del modelo como datos de series temporales separados y modelarlos con otro modelo de pronóstico de series temporales.
Para finalizar el proceso, se utiliza el modelo seleccionado para preveer las futuras observaciones. Como se usaron variables externas como la función tlsm, se debe generar sus valores futuros. Por lo tanto, se sugiere crear un objeto data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras.