#INTRODUCCION
El analisis estadistico es la recompilación e interpretación de datos para cubrir patrones y tendencias. En este informe afianzaremos temas vistos en clase como lo es el modelamiento predictivo de regresión lineal, este tema es interesante ya que es el más utilizado a la hora de predecir valores de una variable cuantitativa apartir de los valores de otra variable cuantitativa. Además, en este informe tambien encontraremos ejercicios de modelamiento de la tendencia de la serie, este modelo es una extensión de los modelos de regresión, aunque es util para describir pautas que siguen una serie temporal. y por ultimo, los componentes estacionales que son fundamentales para el analisis a corto plazo. #Primera parte
Ingeniería de características de los componentes de la serie
Antes de crear las entradas de regresión que representan la tendencia de la serie y componentes estacionales, primero debemos entender su estructura. En los siguientes ejemplos.Demostrar el proceso de creación de nuevas funciones a partir de series existentes utilizando el gas de EE. UU.serie.Primero cargamos la serie desde el paquete TS studio nuevamente y tracemos con ts_plot.
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")Además, recordemos las principales características de la serie aplicando 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 10Como se puede observar en la trama 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. Nosotros podemos explorar la estructura de los componentes de la serie con la función ts_decompose más a fondo:
ts_decompose(USgas)Podemos observar que en el gráfico anterior la tendencia de la serie es bastante plana entre 2000 y 2010, y tiene un crecimiento bastante lineal en el futuro. Por lo 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 lm, la función de regresión lineal R incorporada de las estadísticas paquete, tendremos que transformar la serie de un objeto ts a un objeto data.frame. Por lo tanto, utilizaremos la función ts_to_prophet del paquete TSstudio:
USgas_df <- ts_to_prophet(USgas)La función transforma el objeto ts en dos columnas de data.frame, donde las dos columnas representan el tiempo y los componentes numéricos de la serie, respectivamente:
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 de transformar la serie en un objeto data.frame, se puede empezar a crear las características de entrada de regresión. La primera característica que crearemos es la tendencia de la serie. El enfoque para construir la variable de tendencia es indexar las observaciones de la serie en orden cronológico:
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. Ya que 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, la frecuencia de las unidades representan los meses del año y, por lo tanto, crearemos una variable categórica con 12 categorías, cada categoría correspondiente a un mes específico del año. Usaremos la función de 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, union USgas_df$seasonal <- factor(month(USgas_df$ds, label = T), ordered = FALSE)Usamos la función factorial para convertir la salida de la función mes en no ordenada variable categórica. Revisemos ahora el marco de datos después de agregar las nuevas funciones:
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 retroceder la serie con esas características, dividiremos las series en una partición de entrenamiento y prueba. Estableceremos los últimos 12 meses de la serie como partición de prueba:
h <- 12 # setting a testing partition length
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ] Ahora, después de crear los marcos de datos de entrenamiento y prueba, recordemos cómo la regresión modelo captura cada uno de los componentes por separado y todos juntos.
#Modelando la tendencia de la serie y estacional componentes.
Primero modelamos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, en la partición de entrenamiento:
md_trend <- lm(y ~ trend, data = train)Usamos 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 se puede observar 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 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 anteriormente.
Como se 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 debido que es realmente cierto, o el error de tipo I. Por lo tanto, para el valor p menor que, el valor umbral, rechazaremos la hipótesis nula con un nivel de significancia de, donde los valores típicos de son 0,1, 0,05, 0,01, etc.
Como siempre, es recomendable que coloque 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.
se usará para predecir los valores ajustados y pronosticados del modelo de tendencia que entrenamos antes de:
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)Crearemos una función de utilidad que traza la serie y la salida del modelo, utilizando el paquete plotly:
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 = "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 objeto data.frame que sigue la misma estructura que el uno de los USgas_df (incluida la variable yhat)
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
Establezcamos las entradas de la función plot_lm con la salida del modelo:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")En general, el modelo fue capaz de capturar el movimiento general de la tendencia, pero una línea de la tendencia puede no captar la ruptura estructural de la tendencia que ocurrió alrededor de 2010. Más tarde, en este capítulo, veremos un método avanzado para capturar una tendencia no lineal.
Por último, pero no menos importante, para el análisis de comparación, mediremos la tasa de error del modelo tanto en la formación y los conjuntos de pruebas:
mape_trend <- c (mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.1644088 0.1299951El proceso de modelado y pronóstico del componente estacional sigue el mismo proceso que aplicamos con la tendencia, al hacer una regresión de la serie con la variable estacional que creamos antes de:
md_seasonal <- lm(y ~ seasonal, data = train)Recordemos los detalles del modelo:
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-16Cómo 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 las incluidas en la pendiente valores. Como se puede observar en el resumen de regresión del modelo estacional, todos los coeficientes son estadísticamente significativos. Además, puede observar que el R-cuadrado ajustado del modelo estacional es algo superior con respecto al modelo tendencial (0,78 frente a 0,1).
Antes de trazar el modelo ajustado y los valores de pronóstico con la función plot_lm, actualizaremos los valores de yhat con la función de predicción:
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)Ahora usaremos la función plot_lm para observar el modelo ajustado y los valores de pronóstico:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")Como se puede observar 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, puede observar que la tendencia de la serie ha desaparecido. Antes de agregar los componentes de tendencia y estacional, puntuamos el rendimiento del modelo:
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.08628973 0.21502100La alta tasa de error en el conjunto de prueba está relacionada con el componente de tendencia que no se ha incluido en el modelo. El siguiente paso es unir los dos componentes en un modelo y pronosticar los valores de las características de la serie:
md1 <- lm(y ~ seasonal + trend, data = train)Revisemos el resumen del modelo después de hacer una 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-16La regresión de la serie con los componentes de tendencia y estacional juntos proporcionan elevación adicional a la R-cuadrada ajustada del modelo de 0,78 a 0,91. Esto se puede observar en el gráfico de la salida del modelo:
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend and Seasonal Components of the
Series")Medimos la puntuación MAPE del modelo en las particiones de entrenamiento y prueba:
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.04774945 0.09143290La regresión de la serie con los componentes de tendencia y estacional proporciona una aumento significativo tanto en la calidad de ajuste del modelo como en la precisión del modelo. No obstante, al mirar la gráfica del ajuste y pronóstico del modelo, se logra observar que la tendencia del modelo es demasiado lineal y falta la ruptura estructural de la tendencia de la serie. Este es el punto en el que agregar un componente polinomial para el modelo podría proporcionar una mejora adicional para la precisión del modelo.
Una técnica simple para capturar una tendencia no lineal es agregar un componente polinomial a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Usaremos el argumento I, lo que nos permite aplicar operaciones matemáticas en cualquiera de los objetos de entrada. Por ello, Usaremos 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 manera:
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-16Al agregar el 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 el siguiente gráfico de salida del modelo, este simple cambio en la estructura del modelo nos 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 observar en el modelo que sigue la puntuación de MAPE, la precisión del modelo mejoró significativamente al agregar la tendencia polinomial al modelo de regresión, ya que el error en el conjunto de pruebas se redujo de 9.2% a 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.04212618#La función tslm
Hasta ahora, hemos visto el proceso manual de transformar un objeto en una regresión lineal con el formato del modelo de previsión. La función tslm del paquete de pronóstico nos brinda una función para transformar un objeto ts en un modelo de pronóstico de regresión lineal. Utilizando el tslm, se puede configurar el componente de regresión junto con otras funciones. Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones del gas de EE. UU serie 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 ts_split función:
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testA continuación, se aplicará la misma fórmula que utilizamos para crear la previsión md2 del modelo anterior usando 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))Recordemos ahora md3, la salida de la función tslm, y comparemos 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 notar en el resultado anterior, ambos modelos (md2 y md3) son iguales.
El uso de la función tslm tiene varias ventajas, en comparación con la configuración manual de un modelo de regresión para la serie con la función lm:
Eficiencia: no es necesario 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 ciertos casos, los datos de series de tiempo pueden contener patrones inusuales que se repiten con el tiempo o no. Los siguientes son ejemplos de tales eventos: Valores atípicos: un solo evento o eventos que están fuera de los patrones normales de la serie. Rotura estructural: Un evento 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 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 de 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 ajustar los coeficientes del modelo en consecuencia. Por ejemplo, es posible observar 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 en el año 2010 su crecimiento fue relativamente plano, la pendiente de la tendencia cambió posteriormente, con resultados positivos de crecimiento. En este caso, podemos usar una variable binaria que sea igual a cero para las observaciones 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 observar en el resumen del modelo anterior, la variable de ruptura estructural es estadísticamente significativa, 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. Es necesario tener 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 los valores futuros de esas variables, ya que estos 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 serie con multisecionalidad. En los siguientes ejemplos, utilizaremos la serie UKgrid para demostrar el pronóstico del enfoque de una serie multisecionalidad 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 se encuentra 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, y algunas características principales (por ejemplo, formato de datos, variables, frecuencia, etc.).Gracias a esta función podemos 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. frecuencia 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")Ahora utilizaremos 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 se puede observar, esta serie posee 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.
Ahora utilizaremos la función ts_plot para trazar y revisar la estructura de la serie:
ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")Como se puede observar en el gráfico anterior, la serie tiene una clara tendencia bajista y tiene una cadena de 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.
En el siguiente mapa de calor de la serie desde 2016 se pueden evidenciar esos patrones, utilizando la funcion ts_heatmap del paquete TSstudio:
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")En este mapa de calor se puede observar que la demanda general aumenta durante el invierno. semanas (por ejemplo, semanas del calendario 1 a 12 y semanas 44 a 52). Además, es posible 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 características de la Serie UKdaily
Para capturar los componentes estacionales de la serie, se configura la serie como diaria frecuencia y se crean las siguientes dos características:
Indicador del día de la semana Indicador del mes del año
Además, podemos 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)Debemos recordar 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 generarse a partir de la serie). debido a esto, usaremos las funciones de filtro para eliminar las filas de la tabla en las que falta la variable lag365 (es decir, las primeras 365 observaciones). Revisaremos 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. Usaremos la primera marca de tiempo de la serie, 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) Luego de que transformamos la serie en un objeto ts, podemos retroceder y confirmar la suposición que hicimos sobre el nivel de correlación entre la serie y sus rezagos estacionales con la Función ts_acf (una versión interactiva de la función acf del paquete TSstudio). Revisaremos la correlación de la serie con sus rezagos de los últimos cuatro años:
ts_acf(UK_ts, lag.max = 365 * 4)
El gráfico ACF anterior nos confirma nuestra suposición, y la serie tiene una fuerte relación con los retrasos estacionales, en particular el retraso 365, el primer retraso. Como nota al margen, 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.). No obstante, en general, es recomendable no usar 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. También, puede aumentar el sesgo de pronóstico a medida que utilizamos valores pronosticados como entradas.
Ahora, cuando hemos creado las nuevas funciones para el modelo y hemos establecido 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 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 utilizaremos 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 Igual forma, tenemos que dividir las características que creamos para el modelo de regresión (el características estacionales y de retraso) en una partición de entrenamiento y prueba siguiendo exactamente el mismo orden como usamos para el objeto ts correspondiente. Usaremos el índice data.frame para configurar la funcionalidad de 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
Luego de haber creado las diferentes características para el modelo, ya podemos comenzar la capacitación. proceso del modelo de previsión. Utilizaremos la partición de entrenamiento y comenzaremos a entrenar. siguientes tres modelos:
Modelo de línea de base: regresando la serie con el componente estacional y tendencia utilizando las funciones integradas de la función tslm. mientras 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
Iniciar 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 al modelo no siempre da mejores resultados. Primero comenzaremos por el modelo de linea 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 de la serie, correspondiente a las observaciones de la partición de prueba, utilizando el pronóstico. función:
fc_tslm1 <- forecast(md_tslm1, h = h)Ahora, comparemos el rendimiento del modelo en los conjuntos de entrenamiento y prueba usando el función test_forecast:
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)En la gráfica de desempeño anterior se puede observar 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. 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 se puede observar 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.071899Para mejorar la precisión del modelo intentaremos agregar 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 utilizando 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, usando un pronosticar con el modelo entrenado y evaluar el desempeño del modelo: Luego de observar la gráfica de desempeño anterior del segundo modelo, podemos notar 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,agregaremos la variable de retraso al modelo y repitamos el mismo proceso que antes de:
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 ver en la siguiente gráfica:
test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)Es difícil observar a partir de la gráfica de desempeño del tercer modelo, si hay una diferencia del segundo modelo. Por lo tanto, revisaremos 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 5.07% en el conjunto de prueba.
#Selección de modelo
En este momento seleccionaremos un modelo. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Debido a 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 un 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 continuar 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 igual manera, 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 con solo oportunidad y no porque la variable de rezago 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 en varios períodos de prueba. Dejaré que usted realice pruebas retrospectivas para ambos modelos como un ejercicio divertido.
Para la selección de funciones se pueden aplicar métodos más sólidos, como por pasos,cresta o regresión de lazo.Es recomendable que lea acerca de esos metodos ya que no se encuentran en este libro. En el capítulo 12, Predicción con modelos de aprendizaje automático, exploramos avanzado enfoques de regresión con modelos de aprendizaje automático, que incluyen métodos de selección de características.
En aras de la simplicidad, iremos con 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 el próximo 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, haremos 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-16En el anterior gráfico de resumen de residuos es posible observar que los residuos no son blancos, ya que existe cierta correlació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 pueden explicar la variación de los residuos. El principal desafío con este enfoque es que es difícil identificar las variables externas que pueden explicar la variación de los residuos, y son también factible de pronosticar. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, sin embargo, 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 posible tratar los residuos del modelo como datos de series de tiempo separados y para 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. Cómo 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 las fechas correspondientes de las observaciones pronosticadas, y luego extraer de este objeto las día de la semana y 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, es posible 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")#Resumen
En este capítulo, fueron presentadas 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 Ingeniería de características simples, es posible transformar un problema de pronóstico en una regresión lineal. problema. La principal ventaja del modelo de regresión lineal con respecto a otros modelos tradicionales de series de tiempo es la capacidad del modelo para incorporar variables y factores. Sin embargo, este modelo puede manejar series de tiempo con multisecionalidad. patrones, como vimos con el pronóstico de demanda de electricidad del Reino Unido. Por último, pero no menos importante, Los enfoques de previsión que demostramos en este capítulo serán la base para Modelado con modelos de aprendizaje automático que analizaremos en el Capítulo 12, Previsión con modelos de aprendizaje automático. En el próximo capítulo, presentaremos los métodos de suavizado exponencial, una familia de modelos de previsión basados en un enfoque de media móvil ponderada.
#Segunda parte
df<- read.csv("C:\\Users\\ehern\\Downloads\\database.csv\\database.csv")
data <- 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, frequency = 4)#Modelamiento predictivo con regresión lineal #Se grafica una serie de tiempo asociada a los terremotos en los años desde 1965 a 2016 utilizando la función ts_plot función que ofrece la librería TSstudio.
library(TSstudio)
data(dt)## Warning in data(dt): data set 'dt' not foundts_plot(dt,
title = "Terremotos",
Ytitle = "Magnitud",
Xtitle = "Fecha")#Se revisan ahora las características principales de nuestra serie de tiempo utilizando las función ts_info() de TSstudio.
ts_info(dt)## The dt series is a ts object with 1 variable and 205 observations
## Frequency: 4
## 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-04-01 5.8
## 3 1965-07-01 6.2
## 4 1965-10-01 5.8
## 5 1966-01-01 5.8
## 6 1966-04-01 6.7dt_df$trend <- 1:nrow(dt_df)
library(lubridate)
dt_df$seasonal <- factor(months(dt_df$ds),ordered = FALSE)
head(dt_df)## ds y trend seasonal
## 1 1965-01-01 6.0 1 enero
## 2 1965-04-01 5.8 2 abril
## 3 1965-07-01 6.2 3 julio
## 4 1965-10-01 5.8 4 octubre
## 5 1966-01-01 5.8 5 enero
## 6 1966-04-01 6.7 6 abril#Finalmente, antes de empezar a hacer la regresión de la serie con esas características, dividiremos la serie de tiempo en una partición de entrenamiento y otra de prueba. Estableceremos los últimos 12 meses de la serie como partición de prueba.
h<- 12
train <- dt_df[1:(nrow(dt_df)-h), ]
test <- dt_df [(nrow(dt_df)-h+1):nrow(dt_df),]
tail (train)## ds y trend seasonal
## 188 2011-10-01 6.5 188 octubre
## 189 2012-01-01 5.9 189 enero
## 190 2012-04-01 6.0 190 abril
## 191 2012-07-01 5.8 191 julio
## 192 2012-10-01 5.7 192 octubre
## 193 2013-01-01 6.1 193 eneromd_trend <- lm(y ~ trend, data = train)
summary (md_trend)##
## Call:
## lm(formula = y ~ trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5156 -0.2903 -0.1353 0.1353 2.6862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0214378 0.0683283 88.125 <2e-16 ***
## trend -0.0004507 0.0006108 -0.738 0.462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4728 on 191 degrees of freedom
## Multiple R-squared: 0.002842, Adjusted R-squared: -0.002378
## F-statistic: 0.5444 on 1 and 191 DF, p-value: 0.4615Lo 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(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)
}Para ello, utilizaremos la función summary y revisaremos los detalles del modelo
plot_lm(data = dt_df,
train = train,
test = test,
title = "Terremotos anuales")Se puede observar en el resultado de la regresión que la cuarta columna corresponde al valor de p para cada uno de los coeficientes del modelo. Este valor proporciona la probabilidad de que se rechace la hipótesis nula dado que esta es verdadera, o el error de tipo I. Es decir, para el valor p menor al valor del humbral se rechaza la hipótesis nula con un nivel de significación de α, donde los valores típico son 0.1,0.05,0.01, entre otros. La predicción de los valores ajustados y previstos del modelo de tendencia se realizará mediante md_trend
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend## [1] 0.04996490 0.05499783Luego, se realiza una función de utilidad que traza las series y las salidas del modelo, usando el paquete plotly
md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)##
## Call:
## lm(formula = y ~ seasonal, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5592 -0.2708 -0.1271 0.1408 2.6408
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.92708 0.06768 87.578 <2e-16 ***
## seasonalenero 0.13210 0.09522 1.387 0.167
## seasonaljulio -0.05625 0.09571 -0.588 0.557
## seasonaloctubre 0.12500 0.09571 1.306 0.193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4689 on 189 degrees of freedom
## Multiple R-squared: 0.02948, Adjusted R-squared: 0.01407
## F-statistic: 1.914 on 3 and 189 DF, p-value: 0.1288train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)Cuando nos referimos a argumentos de la función hablamos de
data: Los datos de entrada, un objeto data.frame que sigue la misma estructura que la de df_df (incluyendo la variable yhat) train: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo test: Igualmente, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de predicción title: El título de la parcela, por defecto, establecido como NULL Establezcamos las entradas de la función plot_lm con la salida del modelo
plot_lm(data = dt_df,
train = train,
test = test,
title = "Terremotos anuales")Este modelo fue capaz de captar el movimiento general de la tendencia, aunque una tendencia lineal puede no captar la ruptura estructural de la tendencia que se produjo en torno a 2010. Más adelante, en este capítulo, veremos un método avanzado para captar una tendencia no lineal. Por último, pero no menos importante, para el análisis de comparación, queremos medir la tasa de error del modelo tanto en en los conjuntos de entrenamiento y de prueba. La norma usada se define como
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal## [1] 0.05007408 0.06321639md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)##
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5979 -0.2887 -0.1294 0.1375 2.6040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9713188 0.0893354 66.842 <2e-16 ***
## seasonalenero 0.1325611 0.0953297 1.391 0.166
## seasonaljulio -0.0557892 0.0958198 -0.582 0.561
## seasonaloctubre 0.1259216 0.0958255 1.314 0.190
## trend -0.0004608 0.0006065 -0.760 0.448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4694 on 188 degrees of freedom
## Multiple R-squared: 0.03245, Adjusted R-squared: 0.01186
## F-statistic: 1.576 on 4 and 188 DF, p-value: 0.1823A causa de que se regresó la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías, que son las que están incrustadas con los valores de la pendiente. Como se puede ver en el resumen de la regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. Además, se puede observar que el R-cuadrado ajustado de el 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 de previsión con la función plot_lm, vamos a actualizar los valores de yhat con la función predict
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = dt_df,
train = train,
test = test,
title = "Terremotos anuales")Usando el gráfico anterior, se puede visualizar como el modelo está haciendo un trabajo bastante bueno para capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que la tendencia de la serie no presenta una buen predicción. Antes de añadir tanto la tendencia como los componentes estacionales, puntuaremos el rendimiento del modelo
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1## [1] 0.05023547 0.06273204md2 <- 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.6067 -0.2907 -0.1314 0.1437 2.5966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.986e+00 1.183e-01 50.581 <2e-16 ***
## seasonalenero 1.323e-01 9.559e-02 1.384 0.168
## seasonaljulio -5.579e-02 9.607e-02 -0.581 0.562
## seasonaloctubre 1.259e-01 9.607e-02 1.311 0.192
## trend -9.038e-04 2.445e-03 -0.370 0.712
## I(trend^2) 2.283e-06 1.220e-05 0.187 0.852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4706 on 187 degrees of freedom
## Multiple R-squared: 0.03263, Adjusted R-squared: 0.006764
## F-statistic: 1.262 on 5 and 187 DF, p-value: 0.2824Al presentar este alto índice de error en el conjunto de pruebas está relacionado 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
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = dt_df,
train = train,
test = test,
title = "Terremotos")Cuando hablamos de la regresión de las series con la tendencia y los componentes estacionales se puede observar cómo se proporciona una mejora significativamente tanto la calidad del ajuste del modelo como la precisión del mismo. Sin embargo, al observar el gráfico del ajuste del modelo y la predicción, se puede observar que la tendencia del modelo es demasiado lineal, y no tiene la ruptura estructural de la tendencia de la serie. Este es el punto en el que la adición de un componente polinómico para el modelo podría proporcionar una mejora adicional de la precisión del modelo. Una técnica sencilla para captar una tendencia no lineal consiste en añadir un componente polinómico a la tendencia de la serie para capturar la curvatura de la tendencia en el tiempo. Utilizaremos el argumento I que nos permite aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Por lo tanto, utilizaremos este argumento para añadir un segundo grado del polinomio para la entrada de la tendencia
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2## [1] 0.05029071 0.06296566##La función tslm Mientras que para la adición del polinomio de segundo grado al modelo de regresión no supuso una mejora significativa de la bondad de ajuste del modelo. En el otro modelo, como se puede ver en el siguiente gráfico de salida, este simple cambio en la estructura del modelo nos permite capturar la ruptura estructural de la tendencia en el tiempo
dt_split <- ts_split(dt, sample.out = h)
train.ts <- dt_split$train
test.ts <- dt_split$testlibrary(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.6067 -0.2907 -0.1314 0.1437 2.5966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.118e+00 1.172e-01 52.191 <2e-16 ***
## season2 -1.323e-01 9.559e-02 -1.384 0.1681
## season3 -1.881e-01 9.559e-02 -1.967 0.0506 .
## season4 -6.347e-03 9.559e-02 -0.066 0.9471
## trend -9.038e-04 2.445e-03 -0.370 0.7120
## I(trend^2) 2.283e-06 1.220e-05 0.187 0.8518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4706 on 187 degrees of freedom
## Multiple R-squared: 0.03263, Adjusted R-squared: 0.006764
## F-statistic: 1.262 on 5 and 187 DF, p-value: 0.2824#Modelado de eventos únicos y eventos no estacionales
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] <- 1md3 <- 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.5995 -0.2898 -0.1282 0.1415 2.5918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.069e+00 1.209e-01 50.193 <2e-16 ***
## season2 -8.880e-02 9.473e-02 -0.937 0.3497
## season3 -1.678e-01 9.475e-02 -1.771 0.0781 .
## season4 -1.180e-02 9.472e-02 -0.125 0.9010
## trend 2.456e-03 4.025e-03 0.610 0.5423
## I(trend^2) -9.253e-06 1.608e-05 -0.575 0.5657
## s_break -1.832e-01 1.745e-01 -1.050 0.2951
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4806 on 198 degrees of freedom
## Multiple R-squared: 0.02683, Adjusted R-squared: -0.002665
## F-statistic: 0.9096 on 6 and 198 DF, p-value: 0.489#La serie UKgrid
library(UKgrid)
dts <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
head(dts)## 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 1964584df <- read.csv("C:\\Users\\ehern\\Downloads\\database.csv\\database.csv")
datab <- dplyr::select(df, Depth, Magnitude)
library(zoo)
dts<-ts(data$Magnitude, start = 10, end = 700, freq = 4)
ts_info(dts)## The dts series is a ts object with 1 variable and 2761 observations
## Frequency: 4
## Start time: 10 1
## End time: 700 1Finalmente, 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%
ts_plot(dts ,
title = "Terremotos anuales",
Ytitle = "Magnitud",
Xtitle = "Fecha")ts_heatmap (dt , last = NULL , wday = TRUE , color = "Blues" ,
title = "Earthquakes" , padding = TRUE )library(dplyr)
UKdaily <- dt_df %>%
mutate(wday = wday(ds, label = TRUE),
month = month(ds, label = TRUE),
lag365 = dplyr::lag(y, 365)) %>%
filter(!is.na(lag365)) %>%
arrange(ds)str(UKdaily)## 'data.frame': 0 obs. of 8 variables:
## $ ds : 'Date' num(0)
## $ y : num
## $ trend : int
## $ seasonal: Factor w/ 4 levels "abril","enero",..:
## $ s_break : num
## $ wday : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..:
## $ month : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..:
## $ lag365 : numstart_date <- min(UKdaily$ds)## Warning in min.default(structure(numeric(0), class = "Date"), na.rm = FALSE):
## ningún argumento finito para min; retornando Infstart <- c(year(start_date), yday(start_date))#CONCLUSION A partir del presente informe podemos identificar el modelamiento predictivo con regresión lineal, Modelamiento de la tendencia de la serie y los componentes estacionales. Tanto la parte teórica, como la práctica. Además, del uso del programa Rstudio para la solución de ejercicios haciendo uso de algoritmos, enfocados en la producción anual de gas en Estados Unidos y datos anuales de los terremotos en el mundo. Lo que nos permite relacionar la las prácticas realizadas con nuestro campo laboral. Creando casos hipotéticos que serán útiles para el futuro. Incrementando el desarrollo de nuevas técnicas para el análisis de resultados cuando se tiene una serie de datos geológicos.