PREDICCIONES CON REGRESION LINEAR
Juan Perez,Daniela Albor, Wendy Ortega, Mateo Fontalvo
30/10/2021
INTRODUCCIÓN
En el informe hablaremos de temas como: La función tslm, apreciación de una serie con componentes de multiestacionalidad componentes, Entrenamiento y prueba del modelo de previsión, selección del modelo, Analisis de residuos, Finalización de previsión, entre otros. que sin duda nos enseñar la importancia del análisis estadístico.
PARTE 1
MODELAMIENTO PREDICTIVO CON REGRESIÓN LINEAL
Ingeniería de características de los componentes de la serie
- Utilizamos la librería ‘TSstudio’ específicamente la función ’ts_plot’para comenzar a graficar nuestra serie de tiempo relacionada al consumo de gas en USA
library(TSstudio)## Warning: package 'TSstudio' was built under R version 4.1.1data(USgas)
ts_plot(USgas,
title = "Consumo de gas en USA",
Ytitle = "Billiones pies cubico ",
Xtitle = "Año",
color = "green")- Utilizamos la función ‘ts_info()’ para mirar y revisar las características de nuestra serie de tiempo, esta función pertenece a la librería ‘TSstudio’
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-Visualizamos los primeros datos con la función ‘head’
head(USgas)## [1] 2510.5 2330.7 2050.6 1783.3 1632.9 1513.1- Podemos crear apartir de nuestra serie de tiempo múltiples series de tiempo con función ‘ts_decompose’ de ‘TSstudio’
ts_decompose(USgas)- La secuencia se convierte de un objeto ‘ts’ a un objeto ‘data.frame.’ ya que la función de regresión lineal incorporada en R del paquete ‘stats’, de este modo, utilizaremos la función ‘ts_to_prophet’ del paquete ‘TSstudio’
-Esta función convierte el objeto ts en dos columnas data.frame, donde ambas columnas simbolizan los componentes temporales y numéricos de la secuencia.
USgas_df <- ts_to_prophet(USgas)utilizamos la función ‘head’ para mostrar los primeros 6 datos:
head(USgas_df)## ds y
## 1 2000-01-01 2510.5
## 2 2000-02-01 2330.7
## 3 2000-03-01 2050.6
## 4 2000-04-01 1783.3
## 5 2000-05-01 1632.9
## 6 2000-06-01 1513.1- Ya habiendo cambiado las series en un objeto data.frame, comenzaremos a crear las características de entrada de la regresión. La primera característica que crearemos es la tendencia de la serie, utilizando la función ‘nrow’.
USgas_df$trend <- 1:nrow(USgas_df)- Utilizaremos la función ‘month’ del paquete ‘lubridate’ para extraer el mes del año de la variable de fecha ds:
library(lubridate)## Warning: package 'lubridate' was built under R version 4.1.1##
## 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)- Utilizamos la función ‘head’ para convertir la salida de la función mes en una variable categórica no ordenada.
head(USgas_df)## ds y trend seasonal
## 1 2000-01-01 2510.5 1 ene
## 2 2000-02-01 2330.7 2 feb
## 3 2000-03-01 2050.6 3 mar
## 4 2000-04-01 1783.3 4 abr
## 5 2000-05-01 1632.9 5 may
## 6 2000-06-01 1513.1 6 jun- Por último, estableceremos los últimos 12 meses de la serie como partición de prueba y entrenamiento.
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), ]Modelización de la tendencia de la serie y de los componentes estacionales
-En primer lugar, modelamos la tendencia de la serie mediante la regresión de la serie sobre la variable de tendencia en el segmento de entrenamiento utilizando la función ’ lm ’
md_trend <- lm(y ~ trend, data = train)- Utilizaremos la función de ‘summary’ 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-09-El modelo se utiliza para predecir los valores corregidos en la parte de entrenamiento y los valores predichos en la parte de prueba.
-La función “predict” del paquete “stats”, como su nombre indica, predice los valores de los datos de entrada sobre la base del modelo dado.
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)- A continuación utilizaremos el paquete proporcionado por la librería ‘plotly’ para crear una práctica función para trazar las series y la salida del modelo.
library(plotly)## Warning: package 'plotly' was built under R version 4.1.1## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.1.1##
## 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)
}-Graficamos nuestra predicción de tendencia de la serie, estableciendo las entradas de la función ‘plot_lm’ con la salida del modelo. utilizando la función ’ lm ’
NOTA: compararemos con la grafica siguiente
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicción del componente de tendencia de la serie")- Tasa de error, creamos un vector con la función ‘c()’
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.1644088 0.1299951- Seguimos el mismo procedimiento que modelización y pronosticar los componentes estacionales mediante la regresión de las series sobre las variables estacionales. utilizando la función ’ lm ’
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) 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-16-Actualizaremos nuestros valores de ‘yhat’ con la función ‘predict’, y luego ajustaremos los valores de previsión con la función ‘plot_lm’
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)-Graficamos la Predicción del componente estacional de la serie, con ayuda la función plot_lm para visualizar el modelo ajustado y los valores previstos.
#NOTA: concluimos que este método de predicción es más exacto que el anterior.
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicción del componente estacional de la serie")- Tasa de Error. creamos un vector con la función ‘c()’
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.08628973 0.21502100- Reunimos los modelos de componentes en un modelo y pronóstico característicos de la serie. utilizando la función ’ lm ’
md1 <- lm(y ~ seasonal + trend, data = train)-Revisamos los modelo con la funcion ‘summary’ después de hacer la regresión de la serie con la tendencia y componentes estacionales
summary(md1)##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -514.73 -77.17 -17.70 85.80 336.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 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-16-Graficamos la Predecir la tendencia y los componentes estacionales del Serie, utilizando la función ‘train’ y ‘test’
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predecir la tendencia y los componentes estacionales del
Serie")-Tasa de error. creamos un vector con la función ‘c()’
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.04774945 0.09143290- Creamos un segundo grado del polinomio para la entrada de tendencia con ayuda de la funcion ‘lm()’
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) 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-16- Graficamos la predicción de la tendencia (polinómica) y los componentes estacionales de la serie.
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicción de la tendencia (polinómica) y los componentes estacionales
de la serie")-Tasa de Error. creamos un vector con la funcion ‘c()’
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.04212618La funcion TSLM
Por el momento, hemos desarrollado el proceso manual de transformar un objeto ts a un formato de modelo de predicción de regresión lineal.
- La función tslm del paquete de pronóstico proporciona una función incorporada para transformar un objeto ts en un modelo de pronóstico de regresión lineal junto con otras características.
A continuación repetiremos el ejemplo anterior y se hará la predicción de las últimas 12 observaciones de la serie USgas con la función tslm usando la tendencia, el cuadrado de la tendencia y los componentes estacionales.
-Primero, dividamos la serie en particiones de entrenamiento y prueba usando la función ts_split:
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testProcederemos a aplicar la misma fórmula que usamos para crear el modelo de predicción md2 anterior haciendo uso de la función tslm:
library(forecast)## Warning: package 'forecast' was built under R version 4.1.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoomd3 <- tslm(train.ts ~ season + trend + I(trend^2))- En cuanto a md3, el resultado de la función tslm, se compara 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-16Ambos modelos (md2 y md3) son iguales según el resultado anterior.
ventajas de usar la función tslm, en lugar de configurar manualmente un modelo de regresión para la serie con la función lm:
Eficiencia: no requiere transformar la serie en un objeto data.frame y ingeniería de características
El objeto de salida admite toda la funcionalidad del pronóstico ; funciones de precisión y verificación de residuos y paquetes TSstudio test_forecast y plot_forecast
Modelado de eventos unicos y eventos no estacionales.
- Los datos de las series temporales pueden presentar patrones inusuales,tales como:
- Valores atipicos
- Eventos no estacionales
- Ruptura estructura
- Al no expresarse en el modelo de regresión, este tipo de eventos sesgará los coeficientes -estimados, ya que el modelo ponderará esos tipos de eventos junto con los eventos regulares de la serie. La regresión de un modelo tslm con variables externas requiere un objeto “data.frame” separado con las variables correspondientes. El siguiente ejemplo demuestra el proceso de creación - de una variable binaria externa que es igual a 0 antes del año 2010 y 1 después, usando 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] <- 1- Ahora usaremos la nueva característica para remodelar la serie USgas:
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)- Usemos la función de para summary revisar la salida 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-16Previsión de una serie con componentes de multiestacionalidad componentes - un estudio de caso
El modelo de regresión nos brinda ventajas que no lo hacen los modelos tradicionales de series de tiempo como ARIMA o Holt-Winters, proporcionando gran variedad de opciones de personalizado, y además, moldear y realizar pronósticos acerca de series de tiempo con mayor dificultad como las de multisecionalidad.
A continuación se presentarán diferentes ejemplos en los que se haga uso de la serie UKgrid para realizar la predicción de varias estaciones con un modelo de regresión lineal.
La serie UKgrid
La serie UKgrid representa las necesidades de electricidad de la red nacional del Reino Unido y está disponible en el paquete UKgrid. Esta serie representa una serie temporal de datos de alta frecuencia con una frecuencia de media hora.
- Usamos la función extract_grid del paquete UKgrid para definir la cadena, las características principales (por ejemplo, el formato de los datos, las variables de frecuencia, etc.).
Usaremos la estructura data.frame para llevar la serie a la frecuencia diaria:
library(UKgrid)## Warning: package 'UKgrid' was built under R version 4.1.1UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")Vamos a utilizar la función head para comprobar las variables de la serie:
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 1964584Obtenemos el siguiente resultado:
A continuación se definen las dos variables de la serie: * TIMESTAMP: Un objeto de fecha utilizado como marca de tiempo o índice de la serie
- ND: La demanda neta de electricidad
Usamos la función ts_plot para trazar la estructura de la serie y verificar:
ts_plot(UKdaily,
title = "The UK National Demand fo Electricity",
Ytitle = "MW",
Xtitle = "Year")Obtenemos el siguiente resultado:
Como puede ver en el gráfico anterior, la serie es claramente descendente y tiene un patrón de cadena estacional. Como vimos en el Capítulo 6, Análisis de estacionalidad, esta serie muestra diversos patrones de estacionalidad:
Diario: Un ciclo de 365 días al año
Día de la semana: Un ciclo de 7 días
Mensual: Afecto del clima
La evidencia de estos patrones se puede evidenciar en el siguiente mapa de calor de la serie desde 2016, gracias a la función ts_heatmap de TSstudio (paquete):
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")Obtenemos el siguiente resultado:
Como puede verse en la serie del mapa de calor, la demanda global aumenta durante las semanas de invierno (por ejemplo, la semana 1 a 12 y las semana 44 a 52). Además, el salto de línea se puede observar entre semana, ya que la demanda aumenta entre semana y disminuye los fines de semana.
Preprocesamiento e ingeniería de características de la serie UKdaily
Para captar los componentes estacionales de la serie, definiremos la serie como diaria y crearemos las siguientes dos características:
- Indicador del día de la semana
- Indicador del mes del año
Dado que es razonable asumir (demostraremos este supuesto con la función ACF, luego de convertir la serie en un objeto ts) que la serie tiene una fuerte correlación con los rezagos estacionales, creamos una variable con un retraso de 365 observaciones.
Usaremos el paquete dplyr para crear estas nuevas 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)Recalcando que el costo de usar una variable de retraso con una longitud de N es la pérdida de las primeras N observaciones, usamos las funciones de filtro para anular las filas de la tabla que faltan en la variable las primeras 365 observaciones ( lag365).
Luego de agregar esas nuevas caracteristicas podemos proceder a analizar la estructura de la tabla UKdaily:
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 ...A partir de lo anterior, obtenemos lo siguiente:
Teniendo en cuenta que la entrada de la función tslm debe estar en formato ts, principalmente para la cadena), convertimos la serie en un objeto ts.
Usamos la primera marca de tiempo de la serie y las funciones (el día del año) en lenguaje R, como year y yday del paquete lubridate para determinar el punto de partida del objeto:
start_date <- min(UKdaily$TIMESTAMP)
start <- c(year(start_date), yday(start_date))Anteriormente en el capítulo 3, se menciona el objeto la serie temporal, utilizaremos la función ts del paquete stats 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 retroceder y confirmar la hipótesis que hicimos sobre el nivel de correlación entre la serie y sus retrasos estacionales con la función ts_acf (una versión interactiva de la función acf del paquete TSstudio).
Examinaremos la correlación de la serie con sus retrasos durante los últimos cuatro años:
library(TSstudio)
data(USgas)
acf(UK_ts, lag.max = 365 * 4)
A partir de lo anterior, obtenemos lo siguiente:
El gráfico ACF anterior confirma nuestra hipótesis, y la cadena tiene una fuerte relación con los retrasos por estación, en particular el fue de 365, el primer rezago.
Como nota al margen, tenga la seguridad de que la serie también tiene una fuerte correlación con los retrasos semanales (es decir, retraso 7, 14, 21, etc.). Sin embargo, en general, no se recomienda utilizar retrasos menores que el horizonte de pronóstico (por ejemplo, retraso 7 si el horizonte de pronóstico es 365) porque también debe predecir estos retrasos antes de poder utilizarlos. como entrada para el modelo.
Esto requiere un esfuerzo adicional ya que necesita prever estos retrasos. Además, puede aumentar el sesgo de pronóstico porque usamos valores de pronóstico como entradas.
Ahora que hemos creado las nuevas características para el modelo y configurado el objeto ts, podemos dividir la serie de entrada que creamos y el objeto de característica externa correspondiente “UKdaily” en una partición de entrenamiento y una partición de prueba.
Recordando que nuestro objetivo es predecir los próximos 365 y la longitud de la serie es lo suficientemente grande (2.50 observaciones), podemos proceder a establecer la partición de prueba a la longitud del horizonte de pronóstico: 365 observaciones. Estableceremos h, una variable indicadora, en 365, que será útil para definir las particiones y luego el horizonte de pronóstico:
h <- 365Como hicimos anterioremente, dividiremos la serie en particiones de entrenamiento y prueba con la función ts_split :
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$testDe manera similar, debemos dividir las características que construimos para el modelo de regresión (las características estacionales y de retraso) en dos particiones, una de entrenamiento y otra de prueba, siguiendo exactamente el orden que usamos para el objeto ts correspondiente. Hacemos uso de la funcionalidad de índice data.frame para establecer la tabla UKdaily en las particiones de entrenamiento y prueba:
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]Capacitación y ensayo del modelo de predicción.
Se trabajará con 3 procesos de modelados de predicción.
Modelo de referencia
Modelo multiestacional
Un modelo multiestacional con un retardo estacional Se realizará la comparación de los tres modelos con el uso de la función “test_forecast”
Empezaremos con el modelo de referencia, retrocediendo la serie con su estacionalidad y los componentes de la tendencia.
md_tslm1 <- tslm(train_ts ~ season + trend)- Utilizaremos el modelo entrenado, “md_tslm1”, para pronosticar los próximos 365 días de la serie, correspondiente a las observaciones de la partición de prueba, utilizando la función “forecast”
fc_tslm1 <- forecast(md_tslm1, h = h)- Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y pruebas usando el función test_forecast
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)- En el modelo presentado se puede ver la tendencia de la serie y el dia de la estacionalidad, también según los resultados MAPE como podemos ver en la salida de la siguiente función de precisión, es del 6,09% y7,77% en las particiones de entrenamiento y pruebas, respectivamente.
accuracy(fc_tslm1, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.030002e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set -1.740215e+04 123156.6 96785.27 -1.8735336 7.160573 0.8112374
## ACF1 Theil's U
## Training set 0.5277884 NA
## Test set 0.5106681 1.071899- Ahora intentaremos mejorar la precisión del modelo, añadiendo el día de la semana y el mes de el año características al modelo
md_tslm2 <- tslm(train_ts ~ season + trend + wday + month,
data = train_df)- Ahora vamos a repetir el mismo proceso que antes, mediante el uso de un pronosticar con el modelo entrenado y evaluar el rendimiento del modelo.
fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)test_forecast(actual = UK_ts,
forecast.obj = fc_tslm2,
test = test_ts)accuracy(fc_tslm2, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -1.589161e-12 70032.14 51999.67 -0.1725754 3.158853 0.4358522
## Test set -1.765395e+04 81683.91 65842.88 -1.3694262 4.703788 0.5518836
## ACF1 Theil's U
## Training set 0.7508301 NA
## Test set 0.6120506 0.6898049- vamos a añadir la variable lag al modelo, y repetiremos el proceso anterior.
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365,
data = train_df)fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)- El rendimiento del tercer modelo se puede ver en el siguiente gráfico
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)