Informe #3
Joel Ortegón, Juan Quintero, Juan Gil, Karolina Kecan
1/11/2021
Primera parte
El modelo de regresión linear relaciona una variable dependiente con una o más variables independientes. Entre sus usos se encuentra el análisis predictivo y series de tiempo
La regresión linear simple viene dada por:
\[Y=β_0+β_1X+ε\]
Donde:
\[Y\] variable dependiente
\[β\] estimación del término constante o intercepción
\[β_1\] coeficiente de la variable independiente
\[X\] variable independiente
\[ε\] error aleatorio o ruido
Con la regresión linear es posible:
Identificar estructuras de las series (tendencia, ciclo, estacionalidad, irregularidad), características, patrones y outliners
Transforma los items anteriores en datos de entrada y los entrega con una serie que permita un modelo de predicción
Ingeniería de las componentes de series
Para lograr lo anterior, es necesario primero introducir una serie. Para esto se carga el paquete TSstudio, escoger la serie “USgas” y graficarla.
library(TSstudio)## Warning: package 'TSstudio' was built under R version 4.1.1data("USgas")
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Para ver las características de la serie es usa la función “ts_info”
ts_info(USgas)## The USgas series is a ts object with 1 variable and 238 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2019 10A partir de esto, es notable que es una serie mensual con una línea de tendecia estable y estacionalidad mensual. Para obtener las componentes estructurales de la serie se usa la función ts_descompose
ts_decompose(USgas)La tendencia de la serie es plana entre el 2000 y el 2010, para luego tender un crecimiento un poco linear
Antes de proceder con la regresión linear, es necesario convertir los datos los datos de la serie de ts a un objeto data.frame. Para esto se usa “ts_to_prophet”
USgas_df <- ts_to_prophet(USgas)Ahora es posible generar una tabla con la fecha y el valor númerico (consumo mensual de gas)
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.1Para crear la regresión es necesario crear la variable de tendencia en la tabla anterior en orden cronologíco:
USgas_df$trend <- 1:nrow(USgas_df)Seguidamente, se crea la variable de estacionalidad. Se busca medir la contribución de cada unidad de frecuencia en la oscilación de la serie, por tanto se usa una variable categórica para cada unidad de frencuencia (meses, en este caso). Para extraer los meses de cada año se usa la función “month” del paquete “lubridate”
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)La función “factor” tiene la función de convertir la función “month” en una variable categórica ordenada. Es decir, los meses en formato número a escritos.Se revisa la tabla.
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 junAhora es posible empezar con la regresión. Para ello se deivide en entrenamiento y testeo. Se usan los últomos 12 meses para el testeo.
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), ]Con los data frame creados, se pasará a evidenciar cómo el modelo de regresión procesa cada componente por separa y todos juntos.
Modelamiento de la componentes tendencial y estacional
Primeramente, se modelara la tendencia de la serie mediante la regresión de la serie usando la tendencia como varibale en el entrenamiento. Para esto se usa la función “lm”
md_trend <- lm(y ~ trend, data = train)Para observar los detalles del modelo se usa “summary”
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-09Los coeficientes de la variable de tendencia tiene un nivel de significancia de 0,001. Mientras que el ajuste R-cuadrado es bajo. Lo anterior encuentta sentido debido a que la mayoria de las variaciones de la serie están más relacionadas a el patrón estacional.
Ahora, es posible usar el modelo para predecir los valores que se ajustaron en el entrenamiento (fitted) y pronósticar futuros valores (forecasted) en el testeo. La función “predict” se usa para predecir valores ajustados y pronosticados del modelo que se entrenó.
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)Con el paquete “plotly”, se crean funciones que incluyen la serie y los valores predictivos (output) del modelo.
library(plotly)## Warning: package 'plotly' was built under R version 4.1.1## 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)
}Los argumentos de la función son:
-data: datos de entrada, un objeto “data-fram”
-train: los datos del entrenamiento que se usaron en el modelo (de entrenamiento)
-test: los datos de testeo que fueron usados para evaluar el modelo predictivo
-title: título
Ahora se definen los datos de entrada de la función “plot_lm” con el modelo de resultados (output).
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")En general, el modelo fue capaz de capturar la tendencia, sin embargo, una tendencia lineal puede fallar al momento de captar el cambio estructural del 2010.
Por último, se mide la tasa de error en la fase de entrenamiento y testeo.
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.1299951Ahora, para el proceso de modelamiento y predictividad de la componente de estacional, se aplica el mismo proceso que se utilizó para la tendencia, con la excepción de usar la variable de estacionalidad en la regresión de la serie.
md_seasonal <- lm(y ~ seasonal, data = train)Los detalles del modelo son:
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-16Como se observa, los coeficientes del modelo son estadísticamente significativos. El ajuste R-cuadrado del modelo estacional es mayor respecto al modelo diseñado para la tendencia (0.78 y 0.1).
Antes de fijar el modelo para la estacionalidad con la función “plot_lm”, se modifican los valores de “yhat” con la función “predict”
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)Ahora, se usa “plot_lm” para visualizar el modelo fijado y valores predictivos:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Aunque el modelo captura la estructura de patrón estacional de la serie, la tendencia de la serie no es observable. Antes de añadir las componentes de tendencia y estacinalidad, se analiza el desempeño de modelo.
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.08628973 0.21502100El alto valor de error está relacionado con la componente de tendecia que no fue incluida en el model. Ahora, se unen las dos componentes en un solo modelo.
md1 <- lm(y ~ seasonal + trend, data = train)Se revisa el resumen del modelo luego de regresar la serie con las dos componentes (tendencia y estacionalidad)
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-16De esta manera, se obtiene un ajuste R-cuadrado del modelo de 0.78 a 0.91. Esto se evidencia en el resultado del modelo (output).
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Con estos resultado, es posible medir la puntuación MAPE del modelo, tanto en entrenamiento como en testeo.
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.09143290Regresar las series con las componentes de tendencia y estacionalidad provee mejoras en la acertividad del modelo. Sin embargo, como se mencionó anteriormente, el modelo usado en la tendencia es muy linear, haciendo que se pierdan cambios estructurales. Para mitigar el esto se añade una componente polinomial, la cual permite aumentar la precisión del modelo.
Una técnica simple para capturar una tendencia no lineal es agregar un componente polinómico a la tendencia de la serie para capturar la curvatura de dicha tendencia a lo largo del tiempo. Se utilizará el argumento I, lo que nos permite aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Por lo tanto, se hará uso de este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)El resumen del modelo se puede ver de la siguiente forma:
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-16La adición del polinomio de segundo grado al modelo de regresión no dio lugar a una mejora de la bondad de ajuste del modelo. En el otro modelo, como se puede ver en la siguiente gráfica de salida del modelo, este simple cambio en la estructura del modelo permite capturar la ruptura estructural de la tendencia a lo largo del tiempo:
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the 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 mejora significativamente al agregar la tendencia polinómica al modelo de regresión, como el error en el conjunto de pruebas se redujo del 9,2% al 4,5%:
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2## [1] 0.03688770 0.04212618la función tslm
Hasta ahora, hemos visto el proceso manual de transformar un objeto ts en una regresión lineal. formato de modelo de previsión. La función tslm del paquete de pronóstico proporciona una función incorporada para transformar un objeto ts en un modelo de previsión de regresión lineal. Usando la función tslm, puede establecer el componente de regresión junto con otras características.
Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la USgas. serie con la función tslm utilizando la tendencia, cuadrado de la tendencia, y la estacional Componentes. 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$testA continuación, aplicaremos la misma fórmula que se utilizó para crear la previsión md2 anterior. modelo en el que se utilizó 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))Revisemos ahora md3, la salida de la función tslm, y comparémosla 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 se puede observar en la salida anterior, ambos modelos (md2 y md3) son idénticos.
Hay varias 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 (como el precisión y comprobar funciones residuales) y paquetes TSstudio (como el funciones test_forecast y plot_forecast)
Modelado de eventos únicos y eventos no estacionales
En algunos casos, los datos de series de tiempo pueden contener patrones inusuales que se repiten con el tiempo o no. Los siguientes son ejemplos de dichos eventos:
- Valores atípicos: un solo evento o eventos que están fuera de los patrones normales del serie.
- Rotura estructural: evento significativo que cambia los patrones históricos del 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 el que ocurren cambia de un ciclo a otro. Un ejemplo común de tal evento son las vacaciones de Semana Santa, que ocurren todos los años alrededor de marzo / abril.
Al no expresarse en el modelo de regresión, este tipo de eventos sesgará la estimación coeficientes, ya que el modelo ponderará esos tipos de eventos junto con los eventos regulares de las series. El uso de codificación en caliente, binarias o variables de marca podría ayudar al modelo a ignorar este tipo de eventos o a ajustar los coeficientes del modelo en consecuencia.
Por ejemplo, puede observarse en la gráfica de descomposición de la serie USgas mostrada anteriormente que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Mientras que el crecimiento antes de la El año 2010 fue relativamente plano, la pendiente de la tendencia cambió posteriormente, con resultados positivos crecimiento. En este caso, se puede usar una variable binaria que sea igual a cero para las observaciones realizadas antes del año 2010 y un año después.
La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejemplo 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, utilizando el 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:
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 ver en el resumen del modelo anterior, la variable de ruptura estructural es estadísticamente significativo, con un nivel de 0,01. Asimismo, en el caso de valores atípicos o festivos, La codificación en caliente se puede aplicar estableciendo una variable binaria que sea igual a 1 siempre que un valor atípico o ocurre un evento recurrente no estacional, y 0 en caso contrario.
Tenga en cuenta que, una vez que haya entrenado un modelo de pronóstico con tslm función con el uso de variables externas, tendrá que producir el valores futuros de esas variables, ya que se utilizarán como entrada de la pronóstico.
Pronosticar una serie con multisemporalidad componentes - un estudio de caso
Una de las principales ventajas del modelo de regresión, frente al tradicional tiempo modelos de la serie 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 de tiempo complejas como series con multisecionalidad.
En los siguientes ejemplos, usaremos la serie UKgrid para demostrar el pronóstico enfoque de una serie multisecionalidad con un modelo de regresión lineal.
Las series UKgrid
La serie UKgrid representa la demanda de electricidad de la red nacional en el Reino Unido, y es disponible en el paquete UKgrid. Esta serie representa una serie de datos de alta frecuencia con frecuencia de media hora. Utilizaremos la función extract_grid de UKgrid paquete para definir la serie, características principales (por ejemplo, formato de datos, variables, frecuencia, etc.). Esta función de transformación nos permite agregar la serie frecuencia de media hora a una frecuencia más baja, como por hora, diaria o mensual. Como nuestro objetivo aquí es pronosticar la demanda diaria en los próximos 365 días, configuraremos la serie a diario, utilizando la estructura data.frame:
library(UKgrid)## Warning: package 'UKgrid' was built under R version 4.1.1UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")Usaremos la función head para revisar 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 1964584Como puede ver, esta serie tiene dos variables:
- TIMESTAMP: un objeto de fecha que se utiliza 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:
ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")Como puede ver en el gráfico anterior, la serie tiene una clara tendencia bajista y tiene una cadena patrón estacional. 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: ciclo de 7 días
- Mensual: efecto del clima
La evidencia de esos patrones se puede ver 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 durante las semanas de invierno (por ejemplo, semanas del calendario 1 a 12 y semanas 44 a 52). Además, se peude observar el cambio de la serie durante los días de las semanas, a medida que aumenta la demanda durante los días laborables de la semana y disminuye durante el fin de semana.
Preprocesamiento e ingeniería de funciones de las series UKdaily
Para capturar los componentes estacionales de la serie, configuraremos la serie como diaria frecuencia y cree las siguientes dos características:
- Indicador del día de la semana
- Indicador del mes del año
Además, como es razonable suponer (confirmaremos este supuesto con el ACF una vez que hemos transformado la serie en un objeto ts) que la serie tiene una fuerte correlación con los retrasos estacionales, crearemos una variable de retraso con un retraso de 365 observaciones. Utilizaremos el paquete dplyr para crear esas características:
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 usar una variable de rezago con una longitud de N es la pérdida de las primeras N observaciones (ya que los retrasos de esas observaciones no pueden ser generados a partir de la serie). Por lo tanto, usamos las funciones de filtro para eliminar las filas de la tabla en las que falta la variable lag365 (es decir, las primeras 365 observaciones).
Revisemos 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 a un objeto ts. Se usará la primera marca de tiempo de la serie y el año y yday (el día del año) funciona desde el paquete lubridate para establecer el inicio del objeto punto:
start_date <- min(UKdaily$TIMESTAMP)
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 de las estadísticas paquete para establecer 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 suposición hicimos sobre el nivel de correlación entre la serie y sus rezagos estacionales con la Función ts_acf (una versión interactiva de la función acf del paquete TSstudio). Revisaremos la correlación de la serie con sus rezagos de los últimos cuatro años:
library(TSstudio)
data(UK_ts)## Warning in data(UK_ts): data set 'UK_ts' not foundacf(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 retrasos estacionales, en particular el retraso 365, el primer retraso específicamente.
Como nota a tener en cuenta, puede estar seguro de que la serie también tiene una fuerte correlación con los retrasos semanales (es decir, el retraso 7, 14, 21, etc.). Sin embargo, en general, no se recomienda que utilice retrasos que sean más pequeños que el pronóstico horizonte (por ejemplo, retraso 7, cuando el horizonte de pronóstico es 365), como lo hará tener que pronosticar esos retrasos también para poder utilizarlos como una entrada en el modelo. Esto implica un esfuerzo adicional, ya que tendrá que pronosticar esos retrasos. Además, puede aumentar el sesgo de pronóstico a medida que utilizamos valores pronosticados como entradas.
Ahora, una vez que hemos creado las nuevas funciones para el modelo y hemos establecido el objeto ts, estamos listo 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 los próximos 365 observaciones, y la longitud de la serie es lo suficientemente grande (2540 observaciones), podemos permitirse establecer la partición de prueba a la longitud del horizonte de pronóstico: 365 observaciones. Estableceremos h, una variable indicadora, en 365 y la usaremos para definir las particiones y, más adelante, el horizonte de pronóstico:
h <- 365Como antes, dividiremos la serie en particiones de entrenamiento y prueba con ts_split función:
UKpartitions <- ts_split(UK_ts, sample.out = h)
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 lo mismo orden como usamos para el objeto ts correspondiente. Usaremos el índice data.frame funcionalidad para configurar la tabla UKdaily para entrenar y probar particiones:
train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]Capacitación y prueba del modelo de pronóstico
Después de haber creado las diferentes características para el modelo, estamos listos para comenzar la capacitación. proceso del modelo de previsión. Utilizaremoss la partición de entrenamiento y comenzaremos a entrenar el siguientes tres modelos:
- Modelo de línea de base: regresando la serie con el componente estacional y tendencial utilizando las funciones integradas de la función tslm. A medida que establecemos la frecuencia de la serie en 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
- Modelo multiestacional: sumando el día de la semana y el mes del año indicadores para captar la multisecionalidad de la serie.
- Un modelo multiestacional con un desfase estacional: usar, además del 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 usando la puntuación MAPE
- Visualización de los valores ajustados y pronosticados frente a los valores reales de la serie usando la función test_forecast
Comenzar con un modelo simplista (modelo de línea de base) nos permitirá observar si la adición de nuevas funciones contribuye al rendimiento del modelo o si debemos evitarlo, ya que agregar más características o complejidad a la modelo no siempre da mejores resultados.
Empezaremos por el modelo de línea base, retrocediendo la serie con su estacionalidad y tendencia. componentes:
md_tslm1 <- tslm(train_ts ~ season + trend)A continuación, utilizaremos el modelo entrenado, md_tslm1, para pronosticar los próximos 365 días del serie, correspondiente a las observaciones de la partición de prueba, utilizando el pronóstico función:
fc_tslm1 <- forecast(md_tslm1, h = h)Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y prueba usando la funión test_forecast:
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)Podemos observar en la gráfica de desempeño anterior que el modelo de línea de base está haciendo un gran trabajo capturando tanto la tendencia de la serie como la estacionalidad del día del año. En el otro Por otro lado, no logra capturar la oscilación que se relaciona con el día de la semana. La puntuación MAPE del modelo, como podemos ver en la salida de la siguiente función de precisión, es 6.09% y 7.77% en las particiones de entrenamiento y prueba, respectivamente:
accuracy(fc_tslm1, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -9.030002e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set -1.740215e+04 123156.6 96785.27 -1.8735336 7.160573 0.8112374
## ACF1 Theil's U
## Training set 0.5277884 NA
## Test set 0.5106681 1.071899Ahora intentaremos mejorar la precisión del modelo, agregando el día de la semana y el mes de las características del año para el modelo:
md_tslm2 <- tslm(train_ts ~ season + trend + wday + month,
data = train_df)Como ahora estamos usando características de una fuente de datos externa, tenemos que especificar la entrada datos con el argumento de datos. Ahora repetiremos el mismo proceso que antes, pero usando ahora un pronóstico con el modelo entrenado y evaluar el desempeño del modelo:
Al observar la gráfica de desempeño 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 puede ser observado en el modelo de puntuación MAPE, que se redujo a 2,87% y 5,23% en el entrenamiento y prueba de particiones, respectivamente:
accuracy(fc_tslm2, test_ts)
Por último, pero no menos importante, agregaremos la variable de retraso al modelo y repetiremos el mismo proceso que antes:
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365,
data = train_df)
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)El desempeño del tercer modelo se puede visualizar en la siguiente gráfica:
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)Es complicado observar desde el comportamiento de la gráfica del tercer modelo si hay una diferencia significativa a comparación con el segundo modelo.Por lo que, vamos a revisar la puntuación MAPE del tercer modelo.
accuracy(fc_tslm3, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set -1.175109e-11 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
## ACF1 Theil's U
## Training set 0.7500146 NA
## Test set 0.6094666 0.6925767Los resultados del tercer modelo muestran una pequeña mejora en la precisión del modelo con un 2,81% en el conjunto de entrenamiento y un 5,07% en el conjunto de pruebas.
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. Dado que tanto el segundo como el tercer modelo lograron una puntuación MAPE 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 prueba vale la pena costo de usar la variable de rezago (es decir, la pérdida de 365 observaciones y el costo adicional de una grado de libertad para el modelo).
Eso depende: No hay una respuesta sencilla a esta pregunta. Además, depende del objetivo de la previsión. Se recomienda que considere la siguiente prueba:
- La primera pregunta que debe hacerse en este caso: ¿Es la variable de rezago estadísticamente ¿significativo? Si la variable no es estadísticamente significativa, no tiene sentido continuando la discusión, es mejor descartar la variable. En el caso de la tercer modelo, podemos usar 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 significativo a un nivel de 0,001.
- De manera similar, podemos aplicar una sola prueba ANOVA con la función anova de la paquete de estadísticas y compruebe 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 significativo.
- Backtesting: podría darse el caso de que el tercer modelo sea más preciso 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.
En aras de la simplicidad, usaremos los criterios de precisión y seleccionaremos el tercer modelo para pronosticar la serie. El siguiente paso es reentrenar el modelo en todas las series y pronosticar los próximos 365 días:
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365,
data = UKdaily)Análisis de residuos
Justo antes de finalizar el pronóstico, hagamos un análisis de los residuos del modelo seleccionado usando 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 ver en el gráfico de resumen de residuos anterior, los residuos no son blancos ruido, ya que existe cierta autocorrelación entre las series de residuos y sus rezagos. Este 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 este problema es identificar variables adicionales que puedan explicar la variación en los residuos. El principal desafío con este enfoque es que es difícil identificar variables externas que puedan explicar la variación de los residuales y que también 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 anticipación.
Un enfoque alternativo, cuando los patrones que quedan en los residuos del modelo es tratar los residuos del modelo como datos de series de tiempo separados y modelarlos con otras series de tiempo. modelo de previsión. Un enfoque común es una regresión con error ARIMA, que presentar en el Capítulo 11, Pronóstico con modelos ARIMA.
Finalizando el pronóstico
Finalicemos el proceso y utilicemos el modelo capacitado seleccionado para pronosticar el futuro 365 observaciones. Como usamos variables externas con la función tslm, tendremos que generar sus valores futuros. Esto es relativamente simple, ya que usamos variables deterministas. Por lo tanto, crearemos un objeto data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras. Un enfoque simplista es crear primero el las fechas correspondientes de las observaciones pronosticadas, y luego extraiga de este objeto el día de la semana y el mes del año:
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) +
days(1),
by = "day",
length.out = h))A continuación, podemos utilizar la variable de fecha para crear las variables de día y mes con el paquete lubridate:
UK_fc_df$wday <- factor(wday(UK_fc_df$date, label = TRUE), ordered = FALSE)
UK_fc_df$month <- factor(month(UK_fc_df$date, label = TRUE), ordered =
FALSE)
UK_fc_df$lag365 <- tail(UKdaily$ND, h)Usemos la función de pronóstico para crear el pronóstico:
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 de la Paquete TSstudio:
plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")Parte 2
Lo primero es colocar los datos de la serie de tiempo en un chunk o código r para poder visualizarlos.
df<- read.csv("C:\\Users\\USUARIO\\Downloads\\database.csv")
library(EDAWR)##
## Attaching package: 'EDAWR'## The following object is masked from 'package:dplyr':
##
## stormsdata <- dplyr::select(df, Date, Magnitude)
c <- as.Date(data$Date,format="%d/%m/%Y")
library(zoo)## Warning: package 'zoo' was built under R version 4.1.1##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdt <- ts(data$Magnitude, start = 1965, end = 2016, freq = 12)
library(TSstudio)
data(dt)## Warning in data(dt): data set 'dt' not foundts_plot(dt,
title = "Earthquakes",
Ytitle = "Magnitude",
Xtitle = "Date")Se habla de serie de tiempo df cuando se trata de una serie mensual con un fuerte componente estacional y una línea de tendecia mayormente estable. Para explorar y/o analizar la estructura de cada componente de la serie, se utiliza la función ts_descompose:
ts_info(dt)## The dt series is a ts object with 1 variable and 613 observations
## Frequency: 12
## Start time: 1965 1
## End time: 2016 1ts_decompose(dt)dt_df <- ts_to_prophet(dt)
head(dt_df)## ds y
## 1 1965-01-01 6.0
## 2 1965-02-01 5.8
## 3 1965-03-01 6.2
## 4 1965-04-01 5.8
## 5 1965-05-01 5.8
## 6 1965-06-01 6.7dt_df$trend <- 1:nrow(dt_df)
library(lubridate)
dt_df$seasonal <- factor(months.Date(dt_df$ds), ordered = FALSE)
head(dt_df)## ds y trend seasonal
## 1 1965-01-01 6.0 1 enero
## 2 1965-02-01 5.8 2 febrero
## 3 1965-03-01 6.2 3 marzo
## 4 1965-04-01 5.8 4 abril
## 5 1965-05-01 5.8 5 mayo
## 6 1965-06-01 6.7 6 junioEs posible observar en el gráfico anterior que la tendencia de la serie es bastante variable entre 1965 y 2016. Por tanto, la tendencia general entre ese intervalo de tiempo no es estrecitcamente lineal. Esta idea nos ayudará a definir la entrada de tendencia para el modelo de regresión. Antes de utilizar la función lm de regresión lineal, se debenn transformar las series de un objeto ts a un objeto data.frame. Así que se hará uso de la función ts_to_prohpet del paquete TSstudio.
h <- 12
train <- dt_df[1:(nrow(dt_df) - h), ]
test <- dt_df[(nrow(dt_df) - h + 1):nrow(dt_df),]
md_trend <- lm(y ~ trend, data = train)
summary(md_trend)##
## Call:
## lm(formula = y ~ trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5517 -0.3130 -0.1193 0.1691 2.7055
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.993e+00 3.712e-02 161.430 <2e-16 ***
## trend 9.875e-05 1.068e-04 0.924 0.356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4545 on 599 degrees of freedom
## Multiple R-squared: 0.001424, Adjusted R-squared: -0.0002432
## F-statistic: 0.8541 on 1 and 599 DF, p-value: 0.3558Con esta función se logra transformar el objeto ts en dos columnas de data.frame, las cuales representan tanto las componentes temporales como las númericas, respectivamente, de la serie en cuestión.
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)
library(plotly)
plot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = "Date"),
yaxis = list(title = "Magnitude"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Una vez se haya transformado la serie de tiempo de un objeto tipo data.frame, ya es posible crear características de “entrada” de la regresión. Para este caso, la primera en ser creada será la tendencia de la serie. Un enfoque básico para generar la variable de tendecia es registrar las observaciones de la serie en orden cronológico.
plot_lm(data = dt_df,
train = train,
test = test,
title = "Earthquakes")La siguiente característica a crear es el componente estacional. Al querer medir la contribución de cada unidad de frecuencia a la oscilación de la serie, se utilizará una variable de tio categórica para cada unidad de frecuencia. Para el caso de la serie df, las unidades de frecuencia representan los meses del año, y por lo tanto, se creará una variable categórica con 12 categorías, para lo cual, cada una corresponde a un mes específico del año.
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.05242448 0.06844580Usaremos la función factorial para convertir la salida de la función mes mes en no ordenada variable categórica. Así está el marco de datos después de agregar las nuevas funciones:
md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)##
## Call:
## lm(formula = y ~ seasonal, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.590 -0.316 -0.123 0.160 2.610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.07000 0.06458 93.988 <2e-16 ***
## seasonalagosto -0.05600 0.09133 -0.613 0.540
## seasonaldiciembre -0.05400 0.09133 -0.591 0.555
## seasonalenero -0.09745 0.09088 -1.072 0.284
## seasonalfebrero -0.09160 0.09133 -1.003 0.316
## seasonaljulio -0.12800 0.09133 -1.401 0.162
## seasonaljunio 0.00800 0.09133 0.088 0.930
## seasonalmarzo -0.05200 0.09133 -0.569 0.569
## seasonalmayo 0.02000 0.09133 0.219 0.827
## seasonalnoviembre -0.03000 0.09133 -0.328 0.743
## seasonaloctubre -0.04700 0.09133 -0.515 0.607
## seasonalseptiembre -0.04000 0.09133 -0.438 0.662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4567 on 589 degrees of freedom
## Multiple R-squared: 0.00862, Adjusted R-squared: -0.009895
## F-statistic: 0.4656 on 11 and 589 DF, p-value: 0.9243Para finalizar, antes de empezar a realizar la regresión de la serie con las características antes mencionadas, se dividirá la serie de tiempo en una partición de entrenamiento y otra de prueba. Se establecerán los últimos 12 meses de la serie como partición de prueba.
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)
plot_lm(data = dt_df,
train = train,
test = test,
title = "Earthquakes")Ahora, luego de crear los marcos de datos de entrenamiento y prueba, veremos cómo la regresión modelo captura cada uno de los coponentes por separado y todos juntos. Primero se modelará la tendencia de la serie haciendo una regresión de la serie con la variable de la tendencia, en la partición de entrenamiento:
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.05245147 0.06608669Como puede ver en el resultado de la regresión anterior, el coeficiente de la variable de tendencia es estadísticamente significativo a un nivel de 0,001. Sin embargo, el R-cuadrado ajustado de la regresión es bastante baja, 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 lo vimos en las gráficas anteriormente. Como siempre, se recomienda que se ponga algo de contexto a los números con datos visualización. 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. El predecir función del paquete de estadísticas, como su nombre lo indica, predice los valores de un dato de entrada basado en un modelo dado. Lo usaremos para predecir los valores ajustados y pronosticados del modelo de tendencia que trabajamos aantes.
md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6059 -0.3071 -0.1187 0.1611 2.6379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.041e+00 7.209e-02 83.797 <2e-16 ***
## seasonalagosto -5.640e-02 9.135e-02 -0.617 0.537
## seasonaldiciembre -5.479e-02 9.135e-02 -0.600 0.549
## seasonalenero -9.775e-02 9.090e-02 -1.075 0.283
## seasonalfebrero -9.140e-02 9.135e-02 -1.001 0.317
## seasonaljulio -1.283e-01 9.135e-02 -1.405 0.161
## seasonaljunio 7.802e-03 9.135e-02 0.085 0.932
## seasonalmarzo -5.190e-02 9.135e-02 -0.568 0.570
## seasonalmayo 1.990e-02 9.135e-02 0.218 0.828
## seasonalnoviembre -3.069e-02 9.135e-02 -0.336 0.737
## seasonaloctubre -4.759e-02 9.135e-02 -0.521 0.603
## seasonalseptiembre -4.049e-02 9.135e-02 -0.443 0.658
## trend 9.892e-05 1.074e-04 0.921 0.357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4567 on 588 degrees of freedom
## Multiple R-squared: 0.01005, Adjusted R-squared: -0.01015
## F-statistic: 0.4974 on 12 and 588 DF, p-value: 0.9167Ahora crearemos una función de utilidad que traza la serie y la salida del modelo, utilizando el paquete plotly
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = dt_df,
train = train,
test = test,
title = "Predicting the Trend and Seasonal Components of the
Series")Los argumentos de la función son los siguientes datos: Los datos de entrada, un objeto data.frame que sigue la misma estructura que el uno de los USgas_df (incluida la variable that) train: el conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo prueba: Del mismo modo, el conjunto de pruebas correspondiente que se utilizó para evaluar la modelo de previsión title: el título de la trama, por defecto, establecido en NULL.
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.05228441 0.06808885En general, el modelo fue capaz de capturar el movimiento general de la tendencia, pero una línea la tendencia puede no captar la ruptura estructural de la tendencia que ocurrió alrededor de 2010. En adelante, en este capítulo, veremos un método avanzado para capturar una tendencia no lineal. Por último, para el análisis de comparación, queremos medir la tasa de error del modelo tanto en la formación y los conjuntos de prueba
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
## -0.6106 -0.3074 -0.1171 0.1632 2.6330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.046e+00 8.347e-02 72.439 <2e-16 ***
## seasonalagosto -5.639e-02 9.142e-02 -0.617 0.538
## seasonaldiciembre -5.479e-02 9.143e-02 -0.599 0.549
## seasonalenero -9.786e-02 9.098e-02 -1.076 0.283
## seasonalfebrero -9.140e-02 9.142e-02 -1.000 0.318
## seasonaljulio -1.283e-01 9.142e-02 -1.403 0.161
## seasonaljunio 7.803e-03 9.142e-02 0.085 0.932
## seasonalmarzo -5.190e-02 9.142e-02 -0.568 0.570
## seasonalmayo 1.990e-02 9.142e-02 0.218 0.828
## seasonalnoviembre -3.069e-02 9.142e-02 -0.336 0.737
## seasonaloctubre -4.759e-02 9.142e-02 -0.521 0.603
## seasonalseptiembre -4.049e-02 9.142e-02 -0.443 0.658
## trend 4.121e-05 4.306e-04 0.096 0.924
## I(trend^2) 9.586e-08 6.926e-07 0.138 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4571 on 587 degrees of freedom
## Multiple R-squared: 0.01008, Adjusted R-squared: -0.01184
## F-statistic: 0.4598 on 13 and 587 DF, p-value: 0.9461Lo primero es modelar la tendencia de la serie haciendo uso de la regresión de la serie con la variable de tendencia trend, en la partición de entrenamiento.
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = dt_df,
train = train,
test = test,
title = "Predicting the Trend (Polynomial) and Seasonal Components
of the Series")Finalmente, como podemos ver en el modelo que sigue a la puntuación MAPE, la precisión del modelo mejoró significativamente al añadir la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas se redujo del 9.2% al 4.5%
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2## [1] 0.05230284 0.06849585La función tsml
Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie USgas con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. En primer lugar, vamos a dividir la serie en particiones de entrenamiento y de prueba utilizando la función ts_split …y la función ts_split:
dt_split <- ts_split(dt, sample.out = h)
train.ts <- dt_split$train
test.ts <- dt_split$testA continuación, aplicaremos la misma fórmula que utilizamos para crear el modelo de previsión md2 anterior utilizando la función tslm:
library(forecast)
md3 <- tslm(train.ts ~ season + trend + I(trend^2))
summary(md3)##
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6106 -0.3074 -0.1171 0.1632 2.6330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.948e+00 8.268e-02 71.943 <2e-16 ***
## season2 6.459e-03 9.098e-02 0.071 0.943
## season3 4.596e-02 9.098e-02 0.505 0.614
## season4 9.786e-02 9.098e-02 1.076 0.283
## season5 1.178e-01 9.098e-02 1.294 0.196
## season6 1.057e-01 9.098e-02 1.161 0.246
## season7 -3.043e-02 9.098e-02 -0.335 0.738
## season8 4.147e-02 9.098e-02 0.456 0.649
## season9 5.737e-02 9.098e-02 0.631 0.529
## season10 5.027e-02 9.098e-02 0.553 0.581
## season11 6.717e-02 9.098e-02 0.738 0.461
## season12 4.307e-02 9.098e-02 0.473 0.636
## trend 4.121e-05 4.306e-04 0.096 0.924
## I(trend^2) 9.586e-08 6.926e-07 0.138 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4571 on 587 degrees of freedom
## Multiple R-squared: 0.01008, Adjusted R-squared: -0.01184
## F-statistic: 0.4598 on 13 and 587 DF, p-value: 0.9461Modelado de eventos únicos y no estacionales
La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejemplo muestra el proceso de creación de una variable binaria externa que es igual a 0 antes del año 2010 y a 1 después, utilizando la tabla dt_df:
r <- which(dt_df$ds == as.Date("1965-02-01"))
dt_df$s_break <- ifelse(c(dt_df$ds) >= 2016, 1, 0)
dt_df$s_break[r] <- 1Utilicemos la función de resumen para revisar el resultado del modelo:
library(forecast)
md3 <- tslm(dt ~ season + trend + I(trend^2) + s_break, data = dt_df)
summary(md3)##
## Call:
## tslm(formula = dt ~ season + trend + I(trend^2) + s_break, data = dt_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5943 -0.3069 -0.1195 0.1672 2.6209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.966e+00 8.350e-02 71.451 <2e-16 ***
## season2 -1.901e-02 9.026e-02 -0.211 0.833
## season3 3.556e-02 9.025e-02 0.394 0.694
## season4 7.665e-02 9.026e-02 0.849 0.396
## season5 1.119e-01 9.026e-02 1.239 0.216
## season6 7.650e-02 9.026e-02 0.848 0.397
## season7 -5.299e-02 9.026e-02 -0.587 0.557
## season8 1.346e-02 9.025e-02 0.149 0.881
## season9 2.907e-02 9.025e-02 0.322 0.747
## season10 3.193e-02 9.025e-02 0.354 0.724
## season11 3.872e-02 9.025e-02 0.429 0.668
## season12 1.707e-02 9.025e-02 0.189 0.850
## trend 4.336e-05 7.172e-04 0.060 0.952
## I(trend^2) 5.748e-08 9.667e-07 0.059 0.953
## s_break 6.683e-03 9.347e-02 0.072 0.943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4579 on 598 degrees of freedom
## Multiple R-squared: 0.009974, Adjusted R-squared: -0.0132
## F-statistic: 0.4303 on 14 and 598 DF, p-value: 0.9651Las series UKgrid
Como nuestro objetivo aquí es prever la demanda diaria en los próximos 365 días, pondremos la serie a frecuencia diaria utilizando la estructura data.frame:
df<- read.csv("C:\\Users\\USUARIO\\Downloads\\database.csv")
library(EDAWR)
data <- dplyr::select(df, Date, Magnitude)
c <- as.Date(data$Date,format="%d/%m/%Y")
library(zoo)
dt <- ts(data$Magnitude, start = 1965, end = 2016, freq = 12)
library(TSstudio)
data(dt)## Warning in data(dt): data set 'dt' not foundts_plot(dt,
title = "Earthquakes",
Ytitle = "Magnitude",
Xtitle = "Date")Preprocesamiento e ingeniería de características de la serie UKdaily
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 rezago de 365 observaciones. Utilizaremos el paquete dplyr para crear esas características:
library(dplyr)
UKdaily <- dt_df %>%
mutate(wday = wday(ds, label = TRUE),
month = month(ds, label = TRUE),
lag365 = dplyr::lag(y, 6)) %>%
filter(!is.na(lag365)) %>%
arrange(ds)Revisemos la estructura de la tabla de UKdaily después de añadir estas novedades:
str(UKdaily)## 'data.frame': 607 obs. of 8 variables:
## $ ds : Date, format: "1965-07-01" "1965-08-01" ...
## $ y : num 5.9 6 6 5.8 5.9 8.2 5.5 5.6 6 6.1 ...
## $ trend : int 7 8 9 10 11 12 13 14 15 16 ...
## $ seasonal: Factor w/ 12 levels "abril","agosto",..: 6 2 12 11 10 3 4 5 8 1 ...
## $ s_break : num 0 0 0 0 0 0 0 0 0 0 ...
## $ wday : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 5 1 4 6 2 4 7 3 3 6 ...
## $ month : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 7 8 9 10 11 12 1 2 3 4 ...
## $ lag365 : num 6 5.8 6.2 5.8 5.8 6.7 5.9 6 6 5.8 ...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. Utilizaremos la primera marca de tiempo de la serie y las funciones year y yday (el día del año) del paquete lubridate para establecer el punto de partida del objeto del objeto:
start_date <- min(UKdaily$ds)
start <- c(year(start_date), yday(start_date))Como vimos en el capítulo 3, El objeto serie temporal, utilizaremos la función ts del paquete stats para establecer el objeto ts:
UK_ts <- ts(UKdaily$y, start = start, frequency = 6)After we transform the series into a ts object, we can go back and confirm the assumption we made about the correlation level between the series and its seasonal lags with the ts_acf function (an interactive version of the acf function from the TSstudio package). We will review the correlation of the series with its lags from the past four years:
library(TSstudio)
data(UK_ts)## Warning in data(UK_ts): data set 'UK_ts' not foundacf(UK_ts, lag.max = 6*4)
### Ahora, después de haber creado las nuevas características para el modelo y establecer el objeto ts, estamos estamos listos para dividir la serie de entrada y el objeto de características externas correspondiente que hemos creado (UKdaily) en una partición de entrenamiento y otra de prueba. Como nuestro objetivo es predecir las próximas 365 y la longitud de la serie es lo suficientemente grande (2.540 observaciones), podemos podemos permitirnos fijar la partición de prueba en la longitud del horizonte de previsión: 365 observaciones. Fijaremos h, una variable indicadora, en 365 y la utilizaremos para definir las particiones y, más adelante el horizonte de previsión:
h <- 6We will use the data.frame index functionality to set the UKdaily table to training and testing partitions:
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$test
train_df <- dt_df[1:(nrow(dt_df) -h), ]
test_df <- dt_df[(nrow(dt_df) - h + 1):nrow(dt_df), ]We will start with the baseline model, regressing the series with its seasonal and trend components:
library(forecast)
md_tslm1 <- tslm(train_ts ~ season + trend)
fc_tslm1 <- forecast(md_tslm1,h=h)
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)La puntuación MAPE del modelo, como podemos ver en la salida de la siguiente función de precisión:
accuracy(fc_tslm1, test_ts)## ME RMSE MAE MPE MAPE MASE
## Training set 7.401461e-18 0.4530068 0.3276904 -0.5031674 5.258810 0.721972
## Test set -1.214004e-01 0.5493034 0.4654311 -2.7560871 7.684636 1.025444
## ACF1 Theil's U
## Training set 0.03579129 NA
## Test set -0.10643223 0.8592929ts_heatmap (dt , last = NULL , wday = TRUE , color = "Blues" ,
title = "Earthquakes" , padding = TRUE )Análisis de residuos
Just before we finalize the forecast, let’s do some analysis of the selected model residuals using the checkresiduals function:
library(TSstudio)
data(dt)## Warning in data(dt): data set 'dt' not foundacf(dt, lag.max = 365*51)
checkresiduals(dt)## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.