Informe #3 - Capitulo 9
Raúl A. Castro P.
2024-04-23
#CAPITULO 9
##Previsión 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 variable independiente única (regresión lineal univariada) o múltiple (regresión lineal multivariada). Este modelo tiene una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo y, en particular, la predicción de series temporales.
El objetivo de este capítulo son los métodos y enfoques para pronosticar datos de series de tiempo con regresión lineal. Esto incluye métodos para descomponer y pronosticar los componentes de la serie (por ejemplo, la tendencia y los patrones estacionales), manejar eventos especiales (como valores atípicos y días festivos) y utilizar variables externas como regresores.
Este capitulo cubre los siguientes topicos:
-Enfoques de pronóstico con modelos de regresión lineal.
-Extracción y estimación de los componentes de la serie.
-Manejo de rupturas estructurales, valores atípicos y eventos especiales.
-Series de pronóstico con multiestacionalidad.
##Requerimiento técnico.
Los siguientes paquetes se utilizarán en este capítulo:
-TSstudio: Versión 0.1.4 y superior.
-Plotly: Versión 4.8 y superior.
-dplyr: Versión 0.8.1 y superior.
-lubridate: versión 1.7.4 y superior.
-forecast: versión 8.5 y superiores.
Puedes acceder a los códigos de este capítulo desde el siguiente enlace:
https://github.com/PacktPublishing/Hands-On-Time-Series-Analysis-with-R/tree/master/Chapter09
#LA REGRECIÓN LINEAL
El uso principal del modelo de regresión lineal es cuantificar la relación entre la variable dependiente Y (también conocida como variable de respuesta) y las variables independientes X (también conocidas como variables predictoras, impulsoras o regresadoras) de manera lineal. En otras palabras, el modelo expresa la variable dependiente como una combinación lineal de las variables independientes. Una relación lineal entre las variables dependientes e independientes se puede generalizar mediante las siguientes ecuaciones:
-En el caso de una sola variable independiente, la ecuación es la siguiente:
Y₁ = βo + β₁ * X1,i + Ei
-Para n variables independientes, la ecuación tiene el siguiente aspecto:
Y₁ = βo + β1X1,i + β2X2,i +…+ βnXn,i + Ei
Las variables del modelo para estas ecuaciones son las siguientes:
-i representa el índice de observaciones, i = 1,…, N.
-Yi representa la i observación de la variable dependiente.
-Xj representa el valor i de la variable independiente j, donde j = 1,…, n.
-βo representa el valor del término constante (o intercepto).
-β₁, representa los parámetros correspondientes (o coeficientes) de las j variables independientes.
-Ei, define el término de error, que no es más que toda la información que no fue capturada por las variables independientes para la observación i.
El término lineal, en el contexto de la regresión, se refiere a los coeficientes del modelo, los cuales deben seguir una estructura lineal (ya que esto nos permite construir una combinación lineal a partir de las variables independientes). Por otro lado, las variables independientes pueden seguir una formación tanto lineal como no lineal.
Suponiendo que las ecuaciones anteriores representan la verdadera naturaleza de la relación lineal entre las variables dependientes e independientes, entonces el modelo de regresión lineal proporciona una estimación de esos coeficientes (es decir, βo, β1,…, βn), que puede ser formalizado por las siguientes ecuaciones:
-Para el modelo de regresión lineal univariado, la ecuación es la siguiente:
Y = β + β1Χ1
-Para el modelo de regresión lineal multivariado, la ecuación es la siguiente:
Υ₁ = βο + β1Χ1,i + β2X2,i +…+ BnXn,i
Las variables para estas ecuaciones son las siguientes:
-i representa el índice de observación, i = 1,…, N.
-Yi representa la estimación de la variable dependiente i observación.
-Xji representa el valor i de la variable independiente j, donde j = 1,…, n.
-βo representa la estimación del término constante (o intercepto).
-βι,…, În, son la estimación de los parámetros (o coeficientes) correspondientes de las n variables independientes.
La estimación de los coeficientes del modelo se basa en los dos pasos siguientes:
-Definir una función de costo (también conocida como función de pérdida), estableciendo alguna métrica de error para minimizar.
-Aplicar optimización matemática para minimizar la función de costos.
El enfoque de estimación más común es aplicar los mínimos cuadrados ordinarios (OLS). Método N como técnica de optimización para minimizar la suma de cuadrados de los residuos (i=1). Cuadrar los residuos tiene dos efectos:
-Evita que los valores positivos y negativos se cancelen entre sí al sumarlos.
-El efecto cuadrado proporciona una penalización exponencial para los residuos con mayor distancia, a medida que su costo aumenta.
Existen múltiples técnicas de estimación además del MCO, como la máxima verosimilitud, el método de momentos y el bayesiano. Aunque esos métodos no están dentro del alcance del libro, se recomienda encarecidamente leer y conocer enfoques alternativos.
#ESTIMACIÓN DE COEFICIENTES CON EL MÉTODO MCO
El OLS es un método de optimización simple que se basa en álgebra y cálculo lineal básico, o cálculo matricial. (Esta sección es para conocimientos generales; si no está familiarizado con el cálculo matricial, puede omitir esta sección). El objetivo del MCO es identificar los coeficientes que minimizan la suma de cuadrados de los residuos. Supongamos que el residual de la i observación se define como lo siguiente:
E1 = Y1 - Y1
Luego podemos establecer la función de costo con la siguiente expresión:
ΣE1^2 = (Y1 - Y₁)² + (Y2Y2)²+…+(Yn - Yn)2
Antes de aplicar el método MCO para minimizar la suma de cuadrados de los residuos, por razones de simplicidad, transformaremos el representante de la función de costos en una formación matricial:
Esos conjuntos de matrices representan lo siguiente:
-Vector Y (o matriz N × 1), que representa una variable dependiente con N observaciones.
-X, una matriz de N× n + 1 dimensiones, que representa las n variables independientes correspondientes y un escalar de 1 para el componente de intersección (Bo).
-β, una matriz de dimensiones (n + 1) × 1, que representa los coeficientes del modelo.
-E, una matriz de N × 1 dimensiones, que representa el error correspondiente (o la diferencia entre) el valor real Y y su estimación BX.
El término residual Et no debe confundirse con el término de error Ei. Mientras que el primero representa la diferencia entre la serie Y y su estimación Y, el segundo (término de error) representa la diferencia entre la serie y su valor esperado.
Establezcamos la función de costo usando la forma matricial como la definimos anteriormente:
Ahora comenzaremos a expandir esta expresión usando la fórmula de € como describimos anteriormente:
ecuación = (Y – Xb) (Y – Xb)
A continuación, multiplicaremos los dos componentes (€¹ y €) y abriremos los corchetes:
TYTY-2YT XB + BT X XB
Dado que nuestro objetivo es encontrar el factor que minimice esta ecuación, derivaremos la ecuación con respecto a ẞ y luego la pondremos a cero:
θετε дв (ΥΤΥ-2YT XB + BT XT Xẞ) ав = 0
Resolver esta ecuación producirá el siguiente resultado:
XT XB = XTY
Manipular esta ecuación nos permite extraer la matriz ẞ, la estimación de la matriz de coeficientes B:
β = (XX)-¹XTY
Tenga en cuenta que cambiamos la notación de B β a la salida final, ya que representa la estimación del valor real de B.
Las propiedades clave de la estimación de los coeficientes MCO son las siguientes:
-La característica principal del método de estimación de coeficientes MCO es la insesgación de la estimación de los coeficientes con respecto a los valores reales. En otras palabras, para cualquier dado, ẞ, Ε(β) =
-La recta de regresión muestral siempre pasará por la media de X e Y
-La media de Y, la estimación MCO de la variable dependiente Y, es igual a Y, la media de la variable dependiente
-La media del vector de residuos ê es igual a cero, o norte Σ-1 norte 0
#LOS SUPUESTOS DE MCO
Los principales supuestos del modelo OLS son los siguientes:
-Los coeficientes del modelo deben seguir una estructura lineal (por ejemplo, Y = Bo + B₁ex es un modelo lineal pero Y = X βο + X no lo es).
-No existe una colinealidad perfecta entre las variables independientes X1, X2,…, Xn. En otras palabras, 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 distinta de cero (o no constante).
-El término de error E, condicionado a la matriz de variables independientes X, es una variable independiente e idénticamente distribuida (i.i.d) con media 0 y varianza constante σ”.
-Tanto la variable dependiente como la independiente se extraen de la población en una muestra aleatoria. Esta suposición no se cumple cuando se hacen regresiones de datos de series temporales, ya que normalmente las observaciones tienen cierto grado de correlación. Por lo tanto, este supuesto se relaja cuando se hace una regresión de datos de series de tiempo.
El supuesto de MCO sobre el término de error E-i.i.d. con media 0 y 2 la varianza se viola en muchos casos cuando se trabaja con datos de series de tiempo, ya que, normalmente, existe cierto grado de correlación en los residuos del modelo. Esto es principalmente una indicación de que el modelo de regresión no capturó todos los patrones o información de la serie. En el Capítulo 11, Pronósticos con modelos ARIMA, presentaremos un método para manejar este tipo de casos, modelando el término de error con el modelo ARIMA.
#PREVISIÓN CON REGRESIÓN LINEAL
El modelo de regresión lineal, a diferencia de los modelos tradicionales de series de tiempo como el ARIMA o Holt-Winters, no fue diseñado explícitamente para manejar y pronosticar datos de series de tiempo. Más bien, es un modelo genérico con una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo.
Por tanto, la previsión con un modelo de regresión lineal se basa principalmente en los dos pasos siguientes:
Identificar la estructura de la serie, las características clave, los patrones, los valores atípicos y otras características.
Transformar esas características en variables de entrada y hacer una regresión con la serie para crear un modelo de pronóstico.
Las características principales de un modelo de pronóstico de regresión lineal son los componentes de tendencia y estacional. La siguiente sección se centra en identificar la tendencia de la serie y los componentes estacionales y luego transformarlos en variables de entrada del modelo de regresión.
#PREVISIÓN DE LA TENDENCIA Y LOS COMPONENTES ESTACIONALES
En el Capítulo 5, Descomposición de datos de series temporales, introdujimos los componentes estructurales de la serie: los componentes de tendencia, ciclo, estacional e irregular, y los métodos para descomponerlos con la función de descomposición.
Recuerde que la descomposición de una serie se puede definir mediante una de las siguientes expresiones:
-YT+St+Ct+It, cuando la serie tiene estructura aditiva
-Y = T x S x C x It, cuando la serie tiene estructura multiplicativa X
La explicación es la siguiente:
-Tendencia: Representa el crecimiento de la serie en el tiempo después de ajustar y eliminar los efectos estacionales.
-Estacional: un patrón cíclico recurrente que se deriva directamente de las unidades de frecuencia de la serie (por ejemplo, el mes del año para una serie con una frecuencia mensual).
-Ciclo: un patrón cíclico que no está relacionado con la unidad de frecuencia en serie.
-Irregular: cualquier otro patrón que no sea capturado por los componentes de tendencia, estacional y cíclico.
En aras de la simplicidad, podemos eliminar el componente de ciclo, ya que normalmente se fusiona con el componente de tendencia (o ignoramos el componente de ciclo). Por lo tanto, podemos actualizar y reemplazar las ecuaciones anteriores con las siguientes:
-YT+St+It, cuando la serie tiene estructura aditiva
-Yt=Tt St It, cuando la serie tiene estructura multiplicativa
Ahora podemos transformar esas ecuaciones para un modelo de regresión lineal, modificando la notación de las ecuaciones:
Yt = βο+β₁T + B2 St + E
Where:
-Y represents a time series with n observations
-T, an independent variable with n observations, represents the series trend component
-S, an independent variable with n observations, represents the series seasonal component
-Et, the regression error term, represents the irregular component or any pattern that is not captured by the series trend and seasonal component
-Bo, B1, and B2, represent the model intercept, and coefficients of the trend and seasonal components, respectively
Por conveniencia y en el contexto de trabajar con series de tiempo, cambiaremos la notación de observaciones de i a t, ya que representa la dimensión temporal de la serie.
La transformación de una serie con estructura multiplicativa en una formación de regresión lineal requirió una transformación de la serie en una estructura aditiva. Esto se puede hacer aplicando la transformación logarítmica para ambos lados de las ecuaciones:
log(Y) = log(T₁) + log(St) + log(lt)
Una vez que la serie se transforma en una estructura aditiva, la transformación a una formación de regresión lineal es sencilla y sigue el mismo proceso descrito anteriormente:
log(Y) = β + β₁ log(Tt) + β2 log(St) + €t
#CARACTERISTICAS DE INGENIERÍA DE LOS COMPONENTES DE LA SERIE
Antes de crear los datos de regresión que representan la tendencia de la serie y los componentes estacionales, primero debemos comprender su estructura. En los siguientes ejemplos, haremos demostrar el proceso de creación de nuevas funciones a partir de series existentes utilizando USgas serie.
Carguemos la serie desde el paquete TSstudio nuevamente y grafiquemos con la función ts_plot:
library(TSstudio)## Warning: package 'TSstudio' was built under R version 4.3.3data(USgas)
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Además, repasemos las principales características de la serie usando la función ts_info:
ts_info (USgas)
-The USgas series is a ts object with 1 variable and 227 observations
-Frequency: 12
-Start time: 2000 1
-End time: 2018 11
Como se puede ver en el gráfico de la serie, y como vimos en los capítulos anteriores, USgas es una serie mensual con un fuerte componente estacional mensual y una línea de tendencia bastante estable. Podemos explorar más a fondo la estructura de los componentes de la serie con la función ts_decompose:
-ts_decompose (USgas)
library(TSstudio)
# Load the data
data(USgas)
# Display information about the time series
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# Decompose the time series into its components
ts_decompose(USgas)Puede ver en el gráfico anterior que la tendencia de la serie es bastante plana entre 2000 y 2010, y tiene un crecimiento bastante lineal en el futuro. Por tanto, la tendencia general entre 2000 y 2018 no es estrictamente lineal. Esta es una idea importante que nos ayudará a definir la entrada de tendencia para el modelo de regresión.
Antes de usar la función 1m, la función de regresión lineal R incorporada del paquete de estadísticas, tendremos que transformar la serie de un objeto ts a un dato. objeto de marco. Por lo tanto, utilizaremos la función ts_to_prophet del paquete TSstudio:
library(TSstudio)
# Load the data
data(USgas)
# Transform the time series to a data.frame using ts_to_prophet
USgas_df <- ts_to_prophet(USgas)La función transforma el objeto ts en dos columnas de datos. marco, donde las dos columnas representan los componentes numéricos y de tiempo de la serie. Obtenemos el siguiente resultado:
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.1Después transformamos la serie en datos. objeto de marco, podemos comenzar a crear las características de entrada de regresión. La primera característica que crearemos es la tendencia de la serie. Un enfoque básico para construir la variable de tendencia es indexar las series de observaciones en orden cronológico:
# Agregar la característica de tendencia al objeto data.frame
USgas_df$trend <- 1:nrow(USgas_df)La regresión de la serie con el índice de la serie proporciona una estimación del crecimiento marginal de mes a mes, ya que el índice está en orden cronológico con incrementos constantes.
La segunda característica que queremos crear es el componente estacional. Como queremos medir la contribución de cada unidad de frecuencia a la oscilación de la serie, usaremos una variable categórica para cada unidad de frecuencia. En el caso de la serie USgas, las unidades de frecuencia representan los meses del año, por lo que crearemos una variable categórica con 12 categorías, correspondiendo cada categoría a un mes específico del año. Usaremos la función mes del paquete lubridate para extraer el mes del año de la variable de fecha ds:
library (lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionUSgas_df$seasonal <- factor (month (USgas_df$ds, label = T), ordered = FALSE)Usamos la función de factor para convertir el resultado de la función de mes en una variable categórica no ordenada. Revisemos ahora el marco de datos después de agregar las nuevas características. Obtenemos el siguiente resultado:
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 junPor último, pero no menos importante, antes de comenzar a hacer una regresión de la serie con esas características, dividiremos la serie en una partición de entrenamiento y prueba. Estableceremos los últimos 12 meses de la serie como partición de prueba:
# Set the testing partition length
h <- 12
# Create training and testing partitions
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]Ahora, después de crear los marcos de datos de entrenamiento y prueba, revisemos cómo el modelo de regresión captura cada uno de los componentes por separado y todos juntos.
#MODELADO DE LA TENDENCIA DE LA SERIE Y DE LOS COMPONENTES ESTACIONALES
Primero modelaremos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, en la partición de entrenamiento:
.
# Fit the linear regression model
md_trend <- lm(y ~ trend, data = train)Usaremos la función de resumen para revisar 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-09Como puede verse en el resultado de la regresión anterior, el coeficiente de la variable de tendencia es estadísticamente significativo hasta un nivel de 0,001. Sin embargo, el R cuadrado ajustado de la regresión es bastante bajo, lo que generalmente tiene sentido, ya que la mayor parte de la variación de la serie está relacionada con el patrón estacional, como vimos en los gráficos anteriores.
Como puede observar en el resultado de la regresión anterior, la cuarta columna representa el valor p de cada uno de los coeficientes del modelo. El valor p proporciona la probabilidad de que rechacemos la hipótesis nula dado que en realidad es cierta, o el error tipo I. Por lo tanto, para el valor p menor que a, el valor umbral, rechazaremos la hipótesis nula con un nivel de significancia de a, donde los valores típicos de & son 0,1, 0,05, 0,01, etc.
Como siempre, se recomienda poner algo de contexto a los números con visualización de datos. Por lo tanto, usaremos el modelo que creamos para predecir los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La función de predicción del paquete de estadísticas, como su nombre lo indica, predice los valores de unos datos de entrada en función de un modelo determinado.
Lo usaremos para predecir los valores ajustados y pronosticados del modelo de tendencia que entrenamos antes:
# Predictions on the training partition
train$yhat <- predict(md_trend, newdata = train)
# Predictions on the testing partition
test$yhat <- predict(md_trend, newdata = test)Crearemos una función de utilidad que traza la serie y el resultado del modelo, utilizando el paquete plotly:
library(plotly)## Warning: package 'plotly' was built under R version 4.3.3## 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 = "Year"),
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 dato. objeto frame que sigue la misma estructura que la de USgas_df (incluida la variable yhat)
-Entrenar: el conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo.
-Prueba: Asimismo, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de pronóstico
-Título: el título de la trama, de forma predeterminada, establecido en NULL
Configuremos las entradas de la función plot_1m 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, pero una tendencia lineal puede no capturar la ruptura estructural de la tendencia que ocurrió alrededor de 2010. Más adelante, en este capítulo, veremos un método avanzado para capturar un tendencia no lineal.
Por último, pero no menos importante, para el análisis comparativo, queremos medir la tasa de error del modelo tanto en el conjunto de entrenamiento como en el de prueba. Obtenemos el siguiente resultado:
# Calculate MAPE for the training set
mape_train <- mean(abs(train$y - train$yhat) / train$y) * 100
# Calculate MAPE for the testing set
mape_test <- mean(abs(test$y - test$yhat) / test$y) * 100
# Combine the results into a vector
mape_trend <- c(mape_train, mape_test)
# Display MAPE for both training and testing sets
mape_trend## [1] 16.44088 12.99951El proceso de modelación y pronóstico del componente estacional sigue el mismo proceso que aplicamos con la tendencia, haciendo una regresión de la serie con la variable estacional que creamos antes:
# Fit the linear regression model
md_seasonal <- lm(y ~ seasonal, data = train)Repasemos los detalles del modelo. Obtenemos el siguiente resultado:
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) 2774.28 49.75 55.759 < 2e-16 ***
## seasonalfeb -297.92 70.36 -4.234 3.41e-05 ***
## seasonalmar -479.10 70.36 -6.809 9.77e-11 ***
## seasonalabr -905.28 70.36 -12.866 < 2e-16 ***
## seasonalmay -1088.42 70.36 -15.468 < 2e-16 ***
## seasonaljun -1105.49 70.36 -15.711 < 2e-16 ***
## seasonaljul -939.35 70.36 -13.350 < 2e-16 ***
## seasonalago -914.12 70.36 -12.991 < 2e-16 ***
## seasonalsept -1114.74 70.36 -15.843 < 2e-16 ***
## seasonaloct -1022.21 70.36 -14.527 < 2e-16 ***
## seasonalnov -797.53 71.33 -11.180 < 2e-16 ***
## seasonaldic -256.67 71.33 -3.598 0.000398 ***
## ---
## 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 hacemos una 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, que son aquellas integradas con los valores de pendiente. Como puede ver en el resumen de regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. Además, se puede observar que el R cuadrado ajustado del modelo estacional es algo mayor con respecto al modelo de tendencia (0,78 frente a 0,1).
Antes de trazar el modelo ajustado y los valores pronosticados con la función plot_1m, actualizaremos los valores de yhat con la función de predicción:
# Predictions on the training partition
train$yhat <- predict(md_seasonal, newdata = train)
# Predictions on the testing partition
test$yhat <- predict(md_seasonal, newdata = test)Ahora podemos usar la función plot_lm para visualizar el modelo ajustado y los valores de pronóstico. El resultado se muestra de la siguiente manera:
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 está haciendo un trabajo bastante bueno al capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que falta la tendencia de la serie. Antes de agregar los componentes de tendencia y estacional, calificaremos el desempeño del modelo:
# Calculate MAPE for the training set
mape_train_seasonal <- mean(abs(train$y - train$yhat) / train$y) * 100
# Calculate MAPE for the testing set
mape_test_seasonal <- mean(abs(test$y - test$yhat) / test$y) * 100
# Combine the results into a vector
mape_seasonal <- c(mape_train_seasonal, mape_test_seasonal)
# Display MAPE for both training and testing sets
mape_seasonal## [1] 8.628973 21.502100La alta tasa de error en el conjunto de pruebas está relacionada con el componente de tendencia que no se incluyó en el modelo. El siguiente paso es unir los dos componentes en un modelo y pronosticar los valores característicos de la serie:
# Fit the linear regression model
md1 <- lm(y ~ seasonal + trend, data = train)Repasemos el resumen del modelo después de hacer una regresión de la serie con los componentes tendencial y estacional. Obtenemos el siguiente resultado:
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) 2488.8994 32.6011 76.344 < 2e-16 ***
## seasonalfeb -300.5392 41.4864 -7.244 7.84e-12 ***
## seasonalmar -484.3363 41.4870 -11.674 < 2e-16 ***
## seasonalabr -913.1334 41.4880 -22.010 < 2e-16 ***
## seasonalmay -1098.8884 41.4895 -26.486 < 2e-16 ***
## seasonaljun -1118.5855 41.4913 -26.960 < 2e-16 ***
## seasonaljul -955.0563 41.4936 -23.017 < 2e-16 ***
## seasonalago -932.4482 41.4962 -22.471 < 2e-16 ***
## seasonalsept -1135.6874 41.4993 -27.366 < 2e-16 ***
## seasonaloct -1045.7687 41.5028 -25.198 < 2e-16 ***
## seasonalnov -808.0016 42.0617 -19.210 < 2e-16 ***
## seasonaldic -269.7642 42.0635 -6.413 9.05e-10 ***
## 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 los componentes de tendencia y estacional juntos proporciona un impulso adicional al R cuadrado ajustado del modelo de 0,78 a 0,91. Esto se puede ver en el gráfico de salida del modelo. El resultado se muestra de la siguiente manera:
# Make predictions on the training partition
train$yhat <- predict(md1, newdata = train)
# Make predictions on the testing partition
test$yhat <- predict(md1, newdata = test)
# Plot the combined trend and seasonal components
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend and Seasonal Components of the Series")Midamos la puntuación MAPE del modelo tanto en la partición de entrenamiento como en la de prueba. Obtenemos el siguiente resultado:
# Calculate MAPE for the training set
mape_train_md1 <- mean(abs(train$y - train$yhat) / train$y) * 100
# Calculate MAPE for the testing set
mape_test_md1 <- mean(abs(test$y - test$yhat) / test$y) * 100
# Combine the results into a vector
mape_md1 <- c(mape_train_md1, mape_test_md1)
# Display MAPE for both training and testing sets
mape_md1## [1] 4.774945 9.143290Hacer una regresión de la serie con los componentes de tendencia y estacional proporciona un aumento significativo tanto en la calidad del ajuste del modelo como en la precisión del mismo. Sin embargo, al observar el gráfico del ajuste y pronóstico del modelo, se puede observar que la tendencia del modelo es demasiado lineal y le falta la ruptura estructural de la tendencia de la serie. Este es el punto en el que agregar un componente polinomial al modelo podría proporcionar una mejora adicional en la precisión del modelo.
Una técnica sencilla para capturar una tendencia no lineal es agregar un componente polinómico a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Usaremos el argumento I, que nos permite aplicar operaciones matemáticas en cualquiera de los objetos de entrada. Por lo tanto, usaremos este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:
# Fit the linear regression model
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)El resumen del modelo se puede ver de la siguiente manera. Obtenemos el siguiente resultado:
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) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
## seasonalfeb -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
## seasonalmar -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
## seasonalabr -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
## seasonalmay -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
## seasonaljun -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
## seasonaljul -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
## seasonalago -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
## seasonalsept -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
## seasonaloct -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
## seasonalnov -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
## seasonaldic -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-16Agregar el 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 puede ver en el siguiente gráfico de resultados del modelo, este simple cambio en la estructura del modelo nos permite capturar la ruptura estructural de la tendencia a lo largo del tiempo:
# Make predictions on the training partition
train$yhat <- predict(md2, newdata = train)
# Make predictions on the testing partition
test$yhat <- predict(md2, newdata = test)
# Plot the trend (polynomial) and seasonal components
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend (Polynomial) and Seasonal Components of the Series")Como podemos ver en el modelo que sigue la puntuación MAPE, la precisión del modelo mejoró significativamente al agregar la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas cayó del 9,2% al 4,5%:
# Calculate MAPE for the training set
mape_train_md2 <- mean(abs(train$y - train$yhat) / train$y) * 100
# Calculate MAPE for the testing set
mape_test_md2 <- mean(abs(test$y - test$yhat) / test$y) * 100
# Combine the results into a vector
mape_md2 <- c(mape_train_md2, mape_test_md2)
# Display MAPE for both training and testing sets
mape_md2## [1] 3.688770 4.212618#LA FUNCIÓN TSLM
Hasta ahora, hemos visto el proceso manual de transformar un objeto ts a un formato de modelo de pronóstico de regresión lineal. La función ts1m del paquete de pronóstico proporciona una función incorporada para transformar el objeto ats en un modelo de pronóstico de regresión lineal. Con la función ts1m, puede configurar el componente de regresión junto con otras funciones.
Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie USgas con la función ts1m 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:
# Split the time series data into training and testing partitions
USgas_split <- ts_split(USgas, sample.out = h)
# Access the training and testing partitions
train.ts <- USgas_split$train
test.ts <- USgas_split$testA continuación, aplicaremos la misma fórmula que usamos para crear el modelo de pronóstico md2 anterior usando la función ts1m:
library(forecast)## Warning: package 'forecast' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# Fit the linear regression model
md3 <- tslm(train.ts ~ season + trend + I(trend^2))Ahora revisemos md3, la salida de la función ts1m, y compárela con la salida de md2:
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 observar en el resultado anterior, ambos modelos (md2 y md3) son idénticos. Existen varias ventajas al utilizar la función ts1m, en lugar de configurar manualmente un modelo de regresión para la serie con la función 1m:
-Eficiencia: no requiere transformar la serie a un dato. ingeniería de características y objetos de marco
-El objeto de salida admite todas las funciones del pronóstico (como las funciones de precisión y verificación de residuos) y los paquetes de TSstudio (como las funciones test_forecast y plot_forecast)
#MODELADO DE EVENTOS ÚNICOS Y EVENTOS NO ESTACIONALES
En algunos casos, los datos de series temporales pueden contener patrones inusuales que se repiten o no con el tiempo. Los siguientes son ejemplos de tales eventos:
-Valores atípicos: uno o varios eventos únicos que se salen de los patrones normales de la serie.
-Ruptura estructural: Un hecho 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 este tipo de eventos son las vacaciones de Semana Santa, que ocurren cada año alrededor de marzo/abril.
Al no expresarse en el modelo de regresión, este tipo de eventos sesgará los coeficientes estimados, ya que el modelo ponderará ese tipo de eventos junto con los eventos regulares de la serie. El uso de variables de codificación activa, binarias o de bandera podría ayudar al modelo a ignorar este tipo de eventos o ajustar los coeficientes del modelo en consecuencia.
Por ejemplo, se puede observar en el gráfico descompuesto de la serie de gas estadounidense mostrada anteriormente que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Si bien el crecimiento antes del año 2010 fue relativamente plano, la pendiente de la tendencia cambió después, con un crecimiento positivo. . En este caso, podemos utilizar una variable binaria que sea igual a cero para las observaciones anteriores al año 2010 y un año después.
La regresión de un modelo ts1m con variables externas requiere datos separados. objeto de marco con las variables correspondientes. El siguiente ejemplo demuestra el proceso de creación de una variable binaria externa que es igual antes del año 2010 y 1 después, utilizando la tabla USgas_df:
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] <- 1Ahora usaremos la nueva función para remodelar la serie USgas:
# Fit the linear regression model
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)Usemos la función de resumen para revisar el resultado del modelo:
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-16Como se puede observar en el resumen del modelo anterior, la variable ruptura 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, se puede aplicar codificación activa estableciendo una variable binaria que sea igual a 1 siempre que ocurra un evento atípico o recurrente no estacional, y O en caso contrario.
Tenga en cuenta que, una vez que haya entrenado un modelo de pronóstico con la función tslm con el uso de variables externas, tendrá que producir los valores futuros de esas variables, ya que se utilizarán como entrada del pronóstico.
#PRONÓSTICO DE UNA SERIE CON COMPONENTES MULTIESTACIONALES: UN ESTUDIO DE CASO
Una de las principales ventajas del modelo de regresión, frente a los modelos tradicionales de series temporales como ARIMA o Holt-Winters, es que proporciona una amplia gama de opciones de personalización y nos permite modelar y pronosticar datos de series temporales complejas como series con multiestacionalidad.
En los siguientes ejemplos, utilizaremos la serie UKgrid para demostrar el enfoque de pronóstico de una serie multiestacional con un modelo de regresión lineal.
#LA SERIE UKGRID
La serie UKgrid representa la demanda de electricidad de la red nacional en el Reino Unido y está disponible en el paquete UKgrid. Esta serie representa datos de series temporales de alta frecuencia con frecuencia de media hora. Utilizaremos la función extract_grid del paquete UKgrid para definir la serie, las características principales (por ejemplo, formato de datos, variables, frecuencia, etc.). Esta función de transformación nos permite agregar la frecuencia de la serie desde cada media hora a una frecuencia más baja, como por hora, por día o por mes. Como nuestro objetivo aquí es pronosticar la demanda diaria en los próximos 365 días, estableceremos la serie en frecuencia diaria utilizando los datos. estructura del marco:
# Load the required package
library(UKgrid)## Warning: package 'UKgrid' was built under R version 4.3.3# Extract the UK daily electricity demand data
UKdaily <- extract_grid(type = "data.frame", columns = "ND", aggregate = "daily")Usaremos la función head para revisar las variables de la serie. Obtenemos el siguiente resultado:
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 1964584Como puedes ver, esta serie tiene dos variables:
• TIMESTAMP: un objeto de fecha utilizado como marca de tiempo o índice de la serie.
• ND: La demanda neta de electricidad.
Usaremos la función ts_plot para trazar y revisar la estructura de la serie:
# Load the required package
library(TSstudio)
# Plot the time series
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 bajista y un patrón estacional de cadena. Como vimos en el Capítulo 6, Análisis de estacionalidad, esta serie tiene múltiples patrones de estacionalidad:
-Diariamente: Un ciclo de 365 días al año.
-Día de la semana: un ciclo de 7 días
-Mensual: Efectuado por el clima
Se puede ver evidencia de esos patrones en el siguiente mapa de calor de la serie desde 2016 usando 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 de calor de la serie, la demanda general aumenta a lo largo de las semanas de invierno (por ejemplo, las semanas 1 a 12 del calendario y las semanas 44 a 52). Además, se puede observar el cambio de la serie durante los días de la semana, ya que la demanda aumenta durante los días laborables de la semana, y disminuye durante el fin de semana.
#PREPROCESAMIENTO E INGENIERÍA DE FUNCIONES DE LA SERIE UKDAILY
Para capturar los componentes estacionales de la serie, estableceremos la serie como frecuencia diaria y crearemos las dos funciones siguientes:
-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 función ACF una vez que hayamos transformado la serie en un objeto ts) que la serie tiene una fuerte correlación con los rezagos estacionales, crearemos una variable de rezago con un retraso de 365 observaciones. Utilizaremos el paquete dplyr para crear esas funciones:
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)Recuerde que el costo de utilizar una variable de rezago con una longitud de N es la pérdida de las primeras N observaciones (ya que los rezagos de esas observaciones no pueden generarse a partir de la serie). Por lo tanto, utilizamos las funciones de filtro para eliminar las filas de la tabla en las que falta la variable lag365 (es decir, las primeras 365 observaciones).
Repasemos la estructura de la tabla UKdaily después de agregar esas nuevas características:
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 (al menos para la serie), convertiremos la serie en un objeto ts. Usaremos la primera marca de tiempo de la serie y las funciones año e ydía (el día del año) del paquete lubridate para establecer el punto de partida del objeto:
# Calculate the start date of the dataset
start_date <- min(UKdaily$TIMESTAMP)
# Extract the year and day of the year from the start date
start <- c(year(start_date), yday(start_date))Como vimos en el Capítulo 3, El objeto de serie temporal, usaremos la función ts del paquete de estadísticas para configurar el objeto ts:
UK_ts <- ts(UKdaily$ND,
start = start,
frequency = 365)Después de transformar la serie en un objeto ts, podemos regresar y confirmar la suposición que hicimos sobre el nivel de correlación entre la serie y sus desfases 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 rezagos de los últimos cuatro años:
library(TSstudio)
acf(UK_ts, lag.max = 365 * 4)
El gráfico ACF anterior confirma nuestra suposición y la serie tiene una fuerte relación con los rezagos estacionales, en particular el rezago 365, el primer rezago.
Como nota al margen, puede estar seguro de que la serie también tiene una fuerte correlación con los desfases semanales (es decir, desfases 7, 14, 21, etc.). Sin embargo, en general, no se recomienda utilizar retrasos que sean menores que el horizonte de pronóstico (por ejemplo, retraso 7, cuando el horizonte de pronóstico es 365), ya que también tendrá que pronosticar esos retrasos para poder usarlos. como entrada en el modelo. Esto implica un esfuerzo adicional, ya que habrá que pronosticar esos retrasos. Además, puede aumentar el sesgo del pronóstico ya que utilizamos valores pronosticados como datos de entrada.
Ahora, después de haber creado las nuevas características para el modelo y configurar el objeto ts, estamos listos para dividir la serie de entrada y el objeto de características externas correspondiente que creamos (UKdaily) en una partición de entrenamiento y prueba. Como nuestro objetivo es pronosticar las próximas 365 observaciones y la longitud de la serie es lo suficientemente grande (2540 observaciones), podemos permitirnos establecer la partición de prueba a la longitud del horizonte de pronóstico de 365 observaciones. Estableceremos h, una variable indicadora, en 365 y la usaremos para definir las particiones y, más adelante, el horizonte de previsión:
h <- 365Como antes, dividiremos la serie en particiones de entrenamiento y prueba con la función ts_split:
# Split the time series data into training and testing partitions
UKpartitions <- ts_split(UK_ts, sample.out = h)
# Access the training and testing partitions
train_ts <- UKpartitions$train
test_ts <- UKpartitions$testDe manera similar, tenemos que dividir las características que creamos para el modelo de regresión (las características estacionales y de retraso) en una partición de entrenamiento y prueba siguiendo exactamente el mismo orden que usamos para el objeto ts correspondiente. Usaremos los datos. Funcionalidad de índice de cuadros para configurar la tabla UKdaily para entrenar y probar particiones:
# Create the training data frame
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
# Create the testing data frame
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]#ENTRENAMIENTO Y PRUEBA DEL MODELO DE PRONÓSTICO
Una vez que hayamos creado las diferentes funciones para el modelo, estamos listos para comenzar el proceso de capacitación del modelo de pronóstico. Usaremos 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 ts1m. Como configuramos la frecuencia de la serie 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
-Modelo multiestacional: Agregando 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 del modelo estacional indicadores, 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 prueba utilizando la puntuación MAPE
-Visualizar los valores ajustados y pronosticados versus los valores reales de la serie usando la función test_forecast
Comenzar con un modelo simplista (modelo de referencia) nos permitirá observar si agregar nuevas características contribuye al rendimiento del modelo o si debemos evitarlo, ya que agregar más características o complejidad al modelo no siempre produce mejores resultados.
Comenzaremos con el modelo base, haciendo una regresión de la serie con sus componentes estacionales y de tendencia:
# Fit the time series linear model
md_tslm1 <- tslm(train_ts ~ season + trend)A continuación, utilizaremos el modelo entrenado, md_ts 1m1, para pronosticar los próximos 365 días de la serie, correspondientes a las observaciones de la partición de prueba, usando la función de pronóstico:
# Generate forecasts
fc_tslm1 <- forecast(md_tslm1, h = h)Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y prueba usando el función test_forecast:
# Evaluate the forecast accuracy
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)Podemos observar en el gráfico de rendimiento anterior que el modelo de referencia está haciendo un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año. Por otro lado, no logra capturar la oscilación relacionada con el día de la semana. La puntuación MAPE del modelo, como podemos ver en el resultado de la siguiente función de precisión, es 6,09 % y 7,77 % en las particiones de entrenamiento y prueba, respectivamente. Obtenemos el siguiente resultado:
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.071899Ahora intentaremos mejorar la precisión del modelo agregando las características del día de la semana y el mes del año al modelo:
# Fit the time series linear model
md_tslm2 <- tslm(train_ts ~ season + trend + wday + month, data = train_df)Como ahora utilizamos funciones de una fuente de datos externa, tenemos que especificar los datos de entrada con el argumento de datos. Ahora repetiremos el mismo proceso que antes, utilizando un pronóstico con el modelo entrenado y evaluando el rendimiento del modelo:
Al observar el gráfico de rendimiento anterior del segundo modelo, podemos observar la contribución de las características estacionales en el pronóstico, ya que el segundo modelo fue capaz de capturar tanto la tendencia como los patrones multiestacionales de la serie. Esto también se puede observar en la puntuación del modelo MAPE, que cayó al 2,87% y al 5,23% en las particiones de entrenamiento y prueba, respectivamente:
# Calculate accuracy metrics
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.071899Por último, pero no menos importante, agreguemos la variable de retraso al modelo y repitamos el mismo proceso que antes:
# Fit the time series linear model
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365, data = train_df)
# Generate forecasts
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)El desempeño del tercer modelo se puede ver en el siguiente gráfico. Obtenemos la trama como se ve en el siguiente gráfico:
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 hay una diferencia significativa con respecto al segundo modelo. Por 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 prueba.
#SELECCIÓN DE MODELO
Ahora es el momento de seleccionar un modelo. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Como tanto el segundo como el tercer modelo lograron una puntuación MAPE bastante similar, con una pequeña ventaja para el tercer modelo, deberíamos preguntarnos si una mejora de menos del 0,2% en la tasa de error del conjunto de pruebas vale el costo de usar el variable de rezago (es decir, la pérdida de 365 observaciones y el costo adicional de un grado de libertad para el modelo).
Eso depende.
No hay una respuesta sencilla a esta pregunta. Además, depende del objetivo del pronóstico. Se recomienda que consideres la siguiente prueba:
-La primera pregunta que debe hacerse en este caso: ¿La variable de retraso es estadísticamente ¿significativo? Si la variable no es estadísticamente significativa, no tiene sentido continuar la discusión y es mejor descartarla. 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-05El valor p de la variable lag365 indicó que la variable es estadísticamente significativa a un nivel de 0,001.
-De manera similar, 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 ' ' 1La prueba ANOVA también indicó que la variable lag365 es estadísticamente significativa.
-Backtesting: Podría darse el caso de que el tercer modelo sea más preciso simplemente por casualidad y no porque la variable de retraso contribuya a la precisión del modelo, ya que la diferencia es relativamente pequeña. Por lo tanto, el backtesting de ambos modelos podría ayudar a validar si la contribución de la variable de rezago es consistente durante varios períodos de prueba. Dejaré que usted realice pruebas retrospectivas para ambos modelos como un ejercicio divertido.
Se pueden aplicar métodos más sólidos para la selección de características, como la regresión por pasos, de cresta o de lazo. Aunque esos métodos no están dentro del alcance de este libro, se recomienda leer sobre ellos. En el Capítulo 12, Pronósticos con modelos de aprendizaje automático, exploraremos enfoques de regresión avanzados con modelos de aprendizaje automático, que incluyen métodos de selección de características.
En aras de la simplicidad, seguiremos 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:
# Ajustar el modelo de regresión lineal de series temporales
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365, data = UKdaily)#ANÁLSIS DE RESIDUOS
Justo antes de finalizar el pronóstico, hagamos un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals:
checkresiduals (final_md)
##
## Breusch-Godfrey test for serial correlation of order up to 730
##
## data: Residuals from Linear regression model
## LM test = 3301.1, df = 730, p-value < 2.2e-16Como puede verse en el gráfico de resumen de residuos anterior, los residuos no son ruido blanco, ya que existe cierta autocorrelación entre las series de residuos y sus rezagos. Esto es técnicamente una indicación de que el modelo no capturó todos los patrones o información que existe en la serie. Una forma de abordar esta cuestión es identificar variables adicionales que puedan explicar la variación de los residuos. El principal desafío de este enfoque es que es difícil identificar variables externas que puedan explicar la variación de los residuos y que además sean factibles de pronosticar. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, pero es difícil predecir esos patrones climáticos con un año de antelación.
Un enfoque alternativo, cuando quedan patrones en los residuos del modelo, es tratar los residuos del modelo como datos de series de tiempo separados y modelarlos con otro modelo de pronóstico de series de tiempo. Un enfoque común es una regresión con error ARIMA, que presentaremos en el Capítulo 11, Pronósticos con modelos ARIMA.
#FINALIZANDO EL PRONÓSTICO
Finalicemos el proceso y utilicemos el modelo entrenado seleccionado para pronosticar las 365 observaciones futuras. Como usamos variables externas con la función ts1m, tendremos que generar sus valores futuros. Esto es relativamente simple, ya que utilizamos variables deterministas. Por lo tanto, crearemos un dato. objeto de marco con los valores de wday, mes y lag365 para las próximas 365 observaciones futuras. Un enfoque simplista es crear primero las fechas correspondientes de las observaciones pronosticadas y luego extraer de este objeto el día de la semana y el mes del año:
# Crear un marco de datos para almacenar las fechas de los pronósticos
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + 1,
by = "day",
length.out = h))A continuación, podemos utilizar la variable de fecha para crear las variables día y mes con el paquete lubridate:
# Agregar columna para representar el día de la semana
UK_fc_df$wday <- factor(wday(UK_fc_df$date, label = TRUE), ordered = FALSE)
# Agregar columna para representar el mes
UK_fc_df$month <- factor(month(UK_fc_df$date, label = TRUE), ordered = FALSE)
# Agregar columna para el valor rezagado de 365 días
UK_fc_df$lag365 <- tail(UKdaily$ND, h)Usemos la función de pronóstico para crear el pronóstico:
# Generar pronósticos utilizando el modelo final
UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)Por último, pero no menos importante, trazaremos el pronóstico final con la función plot_forecast del paquete TSstudio. Obtenemos el siguiente resultado:
# Trazar los pronósticos
plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")#RESUMEN
En este capítulo, presentamos las aplicaciones de pronóstico del modelo de regresión lineal. Aunque el modelo de regresión lineal no fue diseñado para manejar datos de series de tiempo, con una ingeniería de características simple podemos transformar un problema de pronóstico 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 la capacidad del modelo para incorporar variables y factores externos. Sin embargo, este modelo puede manejar series temporales con patrones multiestacionales, como vimos con el pronóstico de la demanda de electricidad del Reino Unido. Por último, pero no menos importante, los enfoques de pronóstico que demostramos en este capítulo serán la base para el modelado avanzado con modelos de aprendizaje automático que discutiremos en el Capítulo 12, Pronósticos con modelos de aprendizaje automático.
En el próximo capítulo, presentaremos los métodos de suavizado exponencial, una familia de modelos de pronóstico basados en un enfoque de promedio móvil ponderado.