Introducción

Parte 1

Ingeniería de características de los componentes de la serie

Antes de implementar las entradas de regresión lineal que tienen como función representar la tendencia de una serie, junto a sus componentes estacionales de esta, principalmente se debe comprender la estructura. En los ejemplos implementados se puede comprender el proceso de creación de las características partiendo de series ya existentes y apoyándonos de la serie USgas. La serie será cargada nuevamente a partir del paquete TSstudio y será representada por medio de la función ts_plot:

La salida se muestra de la siguiente manera:

library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.1.1
data(USgas)
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")

También, Se revisará las características principales de la siguiente serie, utilizando nuestra 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 10

En el gráfico de la serie podemos observar y analizar según lo visto en los capítulos pasados, Usgas es una serie que tiene un gran componente estacional mensual y esta presenta una línea de tendencia con una estabilidad bastante alta. También se podrá investigar la estructura que tienen los componentes de la serie, apoyándonos de la función ts_decompose:

La salida se muestra de la siguiente manera:

ts_decompose(USgas)

En el grafico anterior, podemos observar que la serie presenta una tendencia bastante plata en los valores 2000 y 2010, luego presenta un aumento bastante lineal hacia delante. Por lo cual, la tendencia general de esta grafica entre los años 2000 y 2018 no es necesariamente lineal. Es de gran importancia reconocer este dato, ya que nos ayudará a definir la entrada de tendencia que será aplicada para el modelo de regresión. Cabe recordar que antes de utilizar la función lm, función que se basa en la regresión integrada en R del paquete stats, se deberá transformar dicha serie de un objeto ts a on objeto data.frame. De tal manera, la función a utilizar será ts_to_prophet del paquete TSstudio:

USgas_df <- ts_to_prophet(USgas)

La función altera las siglas ts en dos columnas de data.frame, donde estas dos columnas hacen alusión a los componentes tanto numéricos como temporales de la serie.

Se obtiene el siguiente resultado:

head(USgas_df)

Luego de haber transformado las series en un objeto data.frame, se puede comenzar a formas las características de entrada en la regresión. La característica número uno que será creada, es la tendencia de la serie. Que tendrá un enfoque bastante básico para poder armar la variable de tendencia es indexar las observaciones que se tengan de manera cronológica:

USgas_df$trend <- 1:nrow(USgas_df)

La regresión de la serie con el índice que tiene dicha serie, nos proporciona una valoración de crecimiento bastante marginal de mes a mes, esto se debe a que el índice se presenta en un orden cronológico con aumentos de manera constante.

La característica numero dos que se creara es el componente estacional. Con el fin de calcular la contribución de cada unidad de frecuencia e la oscilación de la serie, también se utilizará una variable categórica de cada unidad de frecuencia. Para la serie de USgas, las unidades de frecuencia serán representadas por los meses del año y se creara otra variable categórica con un intervalo de 12 categorías, cada categoría que hace alusión a cada mes del año. Finalmente utilizaremos la función month del paquete libridate para obtener 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)

Para este caso se utiliza la función de tipo factorial para convertir la salida de la función mes en una variable categórica no ordenada. Luego se revisará el marco de datos, después de haber añadido las nuevas funciones.

Se obtiene el siguiente resultado:

head(USgas_df)

Para finalizar, antes de poner en práctica la regresión de la serie con dichas características, se deberá dividir, la numero uno se utilizará para una partición de entrenamiento y la otra es de prueba.

Los 12 últimos meses de la serie serán definidos 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), ]

Después de armar los marcos de datos de entretenimiento y de prueba, lo siguiente es revisar como la función del modelo de regresión va adquiriendo los componentes de manera separada y en grupo.

Modelización de la tendencia de la serie y de los componentes estacionales

Lo principal será modelar la tendencia de la serie practicando una regresión de la serie con la variable de tendencia, encima de la partición de entrenamiento:

md_trend <- lm(y ~ trend, data = train)

Usaremos la función de resumen para analizar los detalles del modelo:

summary(md_trend)
## 
## Call:
## lm(formula = y ~ trend, data = train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -547.2 -307.4 -153.2  333.1 1052.6 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1751.0074    52.6435   33.26  < 2e-16 ***
## trend          2.4489     0.4021    6.09 4.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 394.4 on 224 degrees of freedom
## Multiple R-squared:  0.1421, Adjusted R-squared:  0.1382 
## F-statistic: 37.09 on 1 and 224 DF,  p-value: 4.861e-09

Podemos observar en el resultado de la regresión lineal, el coeficiente que están en la variable tendencia en el área estadística es bastante significativo a un nivel de 0,001. En el R-cuadrado ajustado de la regresión es bajo, lo que en general tiene sentido, debido a que gran parte de la variación de la serie tiene cierta relación con el patrón estacional, como se observa en los gráficos anteriores.

Según lo observado en el resultado de la regresión anterior, la columna cuatro hace alusión al valor p de cada uno de los coeficientes del modelo. El valor p nos indica la probabilidad de que rechacemos la hipótesis nula ya que esta es verdadera, o el error de tipo 1. Para el valor p que es menor que α, luego se rechazara la hipótesis nula con un nivel de significación de α, donde encontramos los valores típicos de α que son 0.1, 0.05, 0.01 etc.

Por lo general, siempre es recomendable contextualizar los números con observación de los datos. Ya que, utilizaremos el modelo que fue creado para pronosticas los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La función predict del paquete stats, como es indicada en su nombre, nos arrija los valores de un dato de entrada apoyándose en un modelo.

Para este caso, se utilizará con el fin de predecir los valores ajustados y previstos del modelo de tendencia que se ha entrenado antes.

train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)

Se creará una función de utilidad con el fin de trazar las series y salidas del modelo, apoyándose del 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':
## 
##     layout
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 = "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:

Las entradas de la función plot_lm será establecidas con la salida del modelo:

La salida se muestra de la siguiente manera:

plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")

En general, el modelo es capaz de percibir el movimiento general de la tendencia, sin embargo, esta tendencia lineal puede que no llegue a captar la ruptura estructural de la tendencia que fue. producto al 2010. Seguidamente, en el siguiente capitulo, se observará un método avanzado para percibir la tendencia no lineal. Para finalizar, y no menos importante, en el análisis de comparación se medirá la tasa de error del modelo, ya sea en los grupos de entrenamiento y en los de prueba:

Se obtiene el siguiente resultado:

mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend
## [1] 0.1644088 0.1299951

Se repasan los detalles del modelo:

md_seasonal <- lm(y ~ seasonal, data = train)

Se obtiene el siguiente resultado:

summary(md_seasonal)
## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -608.98 -162.34  -50.77  148.40  566.89 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2774.28      49.75  55.759  < 2e-16 ***
## seasonalfeb   -297.92      70.36  -4.234 3.41e-05 ***
## seasonalmar   -479.10      70.36  -6.809 9.77e-11 ***
## seasonalabr   -905.28      70.36 -12.866  < 2e-16 ***
## seasonalmay  -1088.42      70.36 -15.468  < 2e-16 ***
## seasonaljun  -1105.49      70.36 -15.711  < 2e-16 ***
## seasonaljul   -939.35      70.36 -13.350  < 2e-16 ***
## seasonalago   -914.12      70.36 -12.991  < 2e-16 ***
## seasonalsept -1114.74      70.36 -15.843  < 2e-16 ***
## seasonaloct  -1022.21      70.36 -14.527  < 2e-16 ***
## seasonalnov   -797.53      71.33 -11.180  < 2e-16 ***
## seasonaldic   -256.67      71.33  -3.598 0.000398 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 216.9 on 214 degrees of freedom
## Multiple R-squared:  0.7521, Adjusted R-squared:  0.7394 
## F-statistic: 59.04 on 11 and 214 DF,  p-value: < 2.2e-16

Haremos una regresión de la variable dependiente con una categórica, este modelo de regresión arma coeficientes para 11 de las 12 categorías, estas son las que se encuentran incrustadas con la pendiente. Como se observa en el resumen de la regresión del modelo estacional, los coeficientes del modelo son estadísticamente valores significativos. También, es posible ver que el R-cuadrado ajustado del modelo estacional es mayor con énfasis al modelo de tendencia (0.78 frente a frente a 0.1)

Antes de armar el modelo ajustado y valores de previsión con la función plot_lm, actualizaremos los valores de yhat con función predict:

train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)

Junto a este podemos utilizar la función plot_lm para analizar el modelo ajustado y los valores pronosticados:

La salida se muestra de la siguiente manera:

plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")

Como se puede ver en el gráfico anterior, 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 falta la tendencia de la serie. Antes de añadir tanto la tendencia como los componentes estacionales, puntuaremos el rendimiento del modelo:

Se obtiene el siguiente resultado:

mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal
## [1] 0.08628973 0.21502100

El elevado porcentaje 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:

md1 <- lm(y ~ seasonal + trend, data = train)

Revisemos el resumen del modelo después de hacer la regresión de la serie con la tendencia y componentes estacionales:

Se obtiene el siguiente resultado:

summary(md1)
## 
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -514.73  -77.17  -17.70   85.80  336.95 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2488.8994    32.6011  76.344  < 2e-16 ***
## seasonalfeb   -300.5392    41.4864  -7.244 7.84e-12 ***
## seasonalmar   -484.3363    41.4870 -11.674  < 2e-16 ***
## seasonalabr   -913.1334    41.4880 -22.010  < 2e-16 ***
## seasonalmay  -1098.8884    41.4895 -26.486  < 2e-16 ***
## seasonaljun  -1118.5855    41.4913 -26.960  < 2e-16 ***
## seasonaljul   -955.0563    41.4936 -23.017  < 2e-16 ***
## seasonalago   -932.4482    41.4962 -22.471  < 2e-16 ***
## seasonalsept -1135.6874    41.4993 -27.366  < 2e-16 ***
## seasonaloct  -1045.7687    41.5028 -25.198  < 2e-16 ***
## seasonalnov   -808.0016    42.0617 -19.210  < 2e-16 ***
## seasonaldic   -269.7642    42.0635  -6.413 9.05e-10 ***
## trend            2.6182     0.1305  20.065  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared:  0.9142, Adjusted R-squared:  0.9094 
## F-statistic: 189.2 on 12 and 213 DF,  p-value: < 2.2e-16

La regresión de la serie de tendencia y los componentes estacionales estos proporcionan un incremento extra del R-cuadrado ajustado del modelo, estos valores van de 0.78 a 0.91. Son posibles obsérvalos en el grafico de los resultados modelo:

La salida se muestra de la siguiente manera:

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")

Se mide las puntuaciones MAPE del modelo en las particiones de entrenamiento y de prueba.

Obtenemos el siguiente resultado:

mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1
## [1] 0.04774945 0.09143290

La regresión de las series con la tendencia y componentes estacionales nos informa una mejora bastante notoria, ya sea en la calidad del ajuste del modelo como es la precisión de este mismo. Sin embargo, al analizar la grafica del ajuste del modelo y previsión, es posible ver que la tendencia de este modelo tiene un comportamiento muy lineal y no presenta ruptura estructural en la tendencia de la serie. En este punto, la adición de un componente polinómico para el modelo proporcionaría una mejora adiciona de la precisión del modelo.

La técnica rápida para obtener una tendencia no lineal, se basa en añadir un componente polinómico a la tendencia de dicha serie con el fin de captar la curvatura de la tendencia que tiene en el tiempo. Nos apoyaremos del argumento, que nos ayuda a aplicar funciones matemáticas sobre cualquiera de los objetos de entrada. De esta manera utilizaremos el argumento para agregar un segundo grado del polinomio para la entrada de la tendencia:

md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)

El resumen del modelo puede verse como sigue:

Obtenemos el siguiente resultado:

summary(md2)
## 
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.635e+03  3.224e+01  81.738  < 2e-16 ***
## seasonalfeb  -3.004e+02  3.540e+01  -8.487 3.69e-15 ***
## seasonalmar  -4.841e+02  3.540e+01 -13.676  < 2e-16 ***
## seasonalabr  -9.128e+02  3.540e+01 -25.787  < 2e-16 ***
## seasonalmay  -1.099e+03  3.540e+01 -31.033  < 2e-16 ***
## seasonaljun  -1.118e+03  3.540e+01 -31.588  < 2e-16 ***
## seasonaljul  -9.547e+02  3.540e+01 -26.968  < 2e-16 ***
## seasonalago  -9.322e+02  3.541e+01 -26.329  < 2e-16 ***
## seasonalsept -1.136e+03  3.541e+01 -32.070  < 2e-16 ***
## seasonaloct  -1.046e+03  3.541e+01 -29.532  < 2e-16 ***
## seasonalnov  -8.001e+02  3.590e+01 -22.286  < 2e-16 ***
## seasonaldic  -2.618e+02  3.590e+01  -7.293 5.95e-12 ***
## trend        -1.270e+00  4.472e-01  -2.840  0.00494 ** 
## I(trend^2)    1.713e-02  1.908e-03   8.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9341 
## F-statistic: 246.1 on 13 and 212 DF,  p-value: < 2.2e-16

La adición de polinomio de grado 2 con respecto al modelo de regresión no nos indica una mejora notoria de la bondad de ajuste del modelo. En siguiente modelo, es posible observar en el siguiente grafico de salida del modelo, este cambio en la estructura del modelo nos otorga la captura de ruptura estructural sobre la tendencia en el 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")

La salida se logra ver de la siguiente manera: Podemos observar en el modelo que sigue a la puntuación MAPE, la precisión de este modelo se mejora significativamente al agregar la tendencia polinómica al modelo de regresión, esto se debe al error en el conjunto de pruebas ya que se redujo del 9.2% al 4.5%:

Obtenemos el siguiente resultado:

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

Por el momento, se logra ver el proceso manual de transformación de un objeto ts a un modelo de previsión de regresión lineal. La función tslm del paquete forecast nos brinda una función integrada para convertir un objeto ts en un modelo de previsión de regresión lineal. Con la función tsml, esta logra establecer el componente de regresión junto con características adicionales.

Se repetirá el ejemplo anterior y predeciremos las 12 ultimas observaciones de la serie USgas con la función tslm con ayuda de la tendencia, el cuadrado de la tendencia y los componentes estacionales. En primer lugar, se dividirá la serie en particiones de entrenamiento y de prueba con ayuda de la función ts_split:

USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$test

Ahora se va aplicar la misma formula que utilizamos con el fin de armar el modelo de previsión md2 utilizando 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 zoo
md3 <- tslm(train.ts ~ season + trend + I(trend^2))

Revisemos ahora md3, la salida de la función tslm, y se compara con la salida de md2:

Obtenemos el siguiente resultado:

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-16

Se logra observar en la salida anterior, los modelos (md2 y md3) con iguales.

Al usar la función tslm se logra observar varias ventajas, en vez de establecer de manera manual un modelo de regresión para las series con la función lm:

Modelización de eventos únicos y eventos no estacionales

En casos específicos, los datos de las series temporales nos ayudan a retener patrones inusuales que se repitan en un lapso del tiempo o no. Los siguientes ejemplos del tipo de eventos:

Al no manifestarse en el modelo de regresión, en estos eventos sesgara los coeficientes estimados, esto se debe a que el modelo ponderara ese tipo de eventos junto con los ventos regulares de la serie. El uso de estas variables de codificación en caliente, binarias o banderas nos podrían ser de gran ayuda al modelo a ignorar estos eventos o también podrían ajustar los coeficientes del modelo en consecuencia.

Un ejemplo de estos es cuando se puede observar en el gráfico de descomposición de la serie USgas mostrando con anterioridad, que la tendencia de la serie presento una ruptura estructural en torno al 2010.

Antes del año 2010 el crecimiento era de manera plana, cuando la pendiente de la tendencia cambio, con un crecimiento positivo. Se pudo utilizar variables de carácter binario igual a cero para las observaciones anteriores a el año 2010 y luego un año después.

La regresuin que es presentada de un modelo tslm con variables de carácter externo que necesitan un objeto data.frame de manera separada con las variables correspondientes. El siguiente ejemplo nos expone el proceso para la creación de una variable binaria externa igual a 0 o antes del 2010 y luego 1 año después, apoyándonos de la tabla USgas_df:

r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse(year(USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1

Ahora utilizaremos la nueva función para remodelar la serie USgas:

md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)

Utilicemos 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-16

Como podemos observar en el resumen del modelo anterior, la variable de ruptura estructural es muy significativo en el área de la estadística, presentando un nivel de 0.01. Para el caso de los valores atípicos o las vacaciones de codificación en caliente se logran aplicar, estableciendo una variable binaria igual a 1 siempre y cuando un valor atípico o un evento no estacional recurrente, y 0 si es de manera contraria.

Se debe tener en cuenta que, luego de haber entrenado el modelo de previsión con la función tslm con la utilización de variables externas, se debe producir los valores futuros de las variables, ya que van a ser utilizadas como una entrada de la previsión.

Previsión de una serie con componentes de multiestacionalidad componentes - un estudio de caso

Las principales ventajas para el modelo de regresión, frente a algunos modelos tradicionales de las series temporales, como son ARIMA o Holt-Winters, es que nos permite una amplia gama de opciones para la personalización y nos ayuda a modelar y pronosticas los datos de series temporales complejas como son las series de multiestacionalidad.

En los ejemplos que serán presentados, se utilizara la serie UKgrid para mostar la previsión de una serie multiestacional con modelos de regresión lineal.

La serie UKgrid

La serie UKgrid hace aluion a la demanda de electricidad de la red nacional del Reino Unido, y se encuetra disponibles en el paquete UKgrid. Esta serie nos muestra una serie de carácter temporal de datos de alta frecuencia con frecuencias que se encuentran entre la media hora. Luego, se procede a utilizar la función extract_grid del paquete UKgrid con el fin de definir la serie, las características principales de estas son: (el formato de los datos, las variables frecuencia, etc.). La función de transformación nos ayuda a agregar la serie de media hora a una frecuencia anterior, igual que la horaria, la diaria o la mensual. El objetivo es prever la demanda diaria en los siguientes 365dias, ya que se pondrá la serie en frecuencia diaria utilizando la siguiente estructura data.frame:

library(UKgrid)
## Warning: package 'UKgrid' was built under R version 4.1.1
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")

Utilizaremos la función head para revisar las variables de la serie:

Obtenemos el siguiente resultado:

head(UKdaily)

Como puede ver, esta serie tiene dos variables:

Utilizaremos la función ts_plot para trazar y revisar la estructura de la serie:

A continuación, se muestra la salida:

ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")

Es posible observar en el grafico anterior, la serie presenta una tendencia clara a la baja y presenta una cadena de patrón estacional. Como se pudo ver en el Cap 6, Análisis estacional, la serie presenta múltiples patrones de estacionalidad:

Le evidencia de estos patrones se logra observar en el siguiente mapa de calor de la serie del año 2016 con ayuda de la función ts_heatmap del paquete TSstudio:

A continuación se muestra la salida:

ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")

Es notorio analizando esto en el mapa de calor de la serie, ya que la demanda global esta en aumento a lo largo de las semanas en el invierno (por ejemplo, en las semanas naturales 1 a 12 y en las 44 a 52). Se logra observar el cambio que tiene la serie durante los días de estas semanas, esto se debe a que la demanda va en aumento durante los días de labor de la semana, y disminuye durante el fin de semana.

Preprocesamiento e ingeniería de características de la serie UKdaily

Para ver los componentes estaciones de la serie, se establecerá la serie como una frecuencia diaria y se crearan las dos características siguientes:

Luego, se puede suponer que (se confirmara esta suposición con ayuda de la función ACF una vez se haya convertido la serie en un objeto ts) ya que la serie tiene una fuerte correlación con los rezagos estacionales, se creara una variable de rezago 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, union
UKdaily <- UKdaily %>%
mutate(wday = wday(TIMESTAMP, label = TRUE),
month = month(TIMESTAMP, label = TRUE),
lag365 = dplyr::lag(ND, 365)) %>%
filter(!is.na(lag365)) %>%
arrange(TIMESTAMP)

Cabe repasar que el coste de utilizar una variable de retardo con una longitud de N es la perdida de las primeras N observaciones (esto se debe a que los rezagos de las observaciones no son posibles de generar a partir de la serie). Es por eso que utilizamos las funciones del filtro para quitar las filas de la tabla en las que falta la variable lag365 (hace referencia a las 365 observaciones).

Revisemos la estructura de la tabla de UKdaily después de añadir estas novedades:

Obtenemos el siguiente resultado:

str(UKdaily)
## 'data.frame':    4939 obs. of  5 variables:
##  $ TIMESTAMP: Date, format: "2006-04-01" "2006-04-02" ...
##  $ ND       : int  1718405 1691341 1960325 2023886 2026204 2008422 1981175 1770440 1749715 2012865 ...
##  $ wday     : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 7 1 2 3 4 5 6 7 1 2 ...
##  $ month    : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ lag365   : int  1920069 1674699 1631352 1916693 1952082 1964584 1990895 2003982 1811436 1684720 ...

Como la entrada de la función tslm debe estar en formato ts (al menos para la serie), convertiremos la serie en un objeto ts. 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$TIMESTAMP)
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$ND,
start = start,
frequency = 365)

Luego de transformar la serie en un objeto ts, es posible volver a confirmar la suposición que hicimos sobre el nivel de correlación entre la serie y sus rezagos estacionales con la función acf. Revisaremos la correlación de la serie con sus rezagos de los últimos cuatro años:

Obtenemos el siguiente resultado:

acf(UK_ts, lag.max = 365 * 4)

De acuerdo con el anterior gráfico ACF se infiere entonces que es correcto suponer sobre el nivel de correlación entre la serie y los rezagos estacionales.

Se puede resaltar de igual manera que la serie guarda correlación con los rezagos semanales. Sin embargo, de manera general, los rezagos más pequeños que el horizonte de previsión pueden ser inútiles dado que requiere trabajo adicional al tener que pronosticar estos rezagos pequeños para usarlos como entrada en el modelo, además de que estos valores pronosticados suponen un aumento en el sesgo de la previsión:

Una vez establecido el objeto ts y finalizada la creación de las nuevas características para el modelo, entonces ya es posible dividir la serie de entrada y el objeto de características externas que corresponde al que ya ha sido creado. Lo siguiente será fijar la partición de prueba en la longitud del horizonte de previsión, esto con el fin de poder predecir las próximas 365, y dado que la longitud de serie es lo suficientemente grande, se puede realizar sin complicaciones. Lo siguiente es definir las particiones con una variable indicadora h, la misma que más adelante se utilizará para definir el horizonte de previsión:

h <- 365

Nuevamente se dividirá la serie en una partición de entrenamiento y otra de prueba usando la función ts_split, utilizada para dividir la serie en estas dos partes:

UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$test

De igual manera, las características creadas se tienen que dividir para el modelo de regresión (los rasgos estacionales y el de retardo), siguiendo exactamente el mismo orden utilizado para el objeto ts correspondiente, en una partición de entrenamiento y otra de prueba. Para establecer la tabla UKdaily en las particiones de entrenamiento y prueba se utiliza la funcionalidad del índice data.frame:

train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]

Entrenamiento y prueba del modelo de previsión

Ya creadas las diferentes características del modelo, entonces el proceso del modelo de previsión puede iniciar. Para este se utilizará la partición de entrenamiento, además de comenzar a entrenar los siguientes modelos:

Los criterios de comparación de estos modelos se basará en los siguientes criterios: - Rendimiento del modelo en el conjunto de entrenamiento y prueba mediante la puntuación MAPE - Visualización de los valores ajustados y previstos frente a los valores reales de la serie mediante la función test_forecast

Empezando con un modelo de referencia es posible observar si añadir nueva características contribuirá al rendimiento del modelo o si se debe evitar, ya que no siempre se tiene mejores resultados al añadirle características o complejidad al modelo.

Comenzaremos haciendo una regresión lineal al modelo de referencia, con sus componentes estacionales y de tendencia:

md_tslm1 <- tslm(train_ts ~ season + trend)

A continuación, utilizaremos el modelo entrenado, md_tslm1, para prever los próximos 365 días de la serie, correspondientes a las observaciones de la partición de prueba, utilizando la función de previsión de previsión:

fc_tslm1 <- forecast(md_tslm1, h = h)

Comparemos el rendimiento del modelo en los conjuntos de entrenamiento y de prueba utilizando la función función test_forecast:

La salida se muestra de la siguiente manera:

test_forecast(actual = UK_ts,
forecast.obj = fc_tslm1,
test = test_ts)

En el gráfico de rendimiento anterior se observa que el modelo de referencia hace un buen trabajo para aprehender la tendencia de la serie y la estacionalidad del día del año. Por otra parte, no alcanza a captar la oscilación relacionada con el día de la semana. Y de acuerdo con la puntuación MAPE del modelo, la salida de la siguiente función de precisión muestra que las particiones de entrenamiento y prueba son de 6,09% y 7,77%, respectivamente:

Se obtiene el siguiente resultado:

accuracy(fc_tslm1, test_ts)
##                         ME     RMSE       MAE        MPE     MAPE      MASE
## Training set -9.030002e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set     -1.740215e+04 123156.6  96785.27 -1.8735336 7.160573 0.8112374
##                   ACF1 Theil's U
## Training set 0.5277884        NA
## Test set     0.5106681  1.071899

Ahora, añadiendo al modelo las características del día de la semana y el mes del año, se intentará mejorar su precisión:

md_tslm2 <- tslm(train_ts ~ season + trend + wday + month,
data = train_df)

Dado que en este punto se están utilizando características de una fuente de datos externa, se debe especificar los datos de entrada usando el argumento data. En seguida, utilizando una con el modelo entrenado y evaluando el rendimiento del modelo, repetiremos este mismo proceso:

De acuerdo al gráfico de rendimiento anterior del segundo modelo, se pueden observar la contribución de las características estacionales a la previsión, puesto que el segundo modelo pudo capturar la tendencia y los patrones multiestacionales de la serie. Es posible verlo también el la puntuación MAPE del modelo, en el que las particiones de entrenamiento y de prueba se redujeron 2,87% y 5,23%, respectivamente:

Obtenemos el siguiente resultado:

accuracy(fc_tslm1, test_ts)
##                         ME     RMSE       MAE        MPE     MAPE      MASE
## Training set -9.030002e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set     -1.740215e+04 123156.6  96785.27 -1.8735336 7.160573 0.8112374
##                   ACF1 Theil's U
## Training set 0.5277884        NA
## Test set     0.5106681  1.071899

Por último, se añade la variable de retardo al modelo, y se repite el mismo proceso anterior:

md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365,
data = train_df)
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)

A continuación, el gráfico de rendimiento del tercer modelo:

Obtenemos la trama como se ve en el siguiente gráfico:

test_forecast(actual = UK_ts,
forecast.obj = fc_tslm3,
test = test_ts)

En la gráfica, observar una diferencia significativa entre el segundo y el tercer modelo puede tornarse difícil, es por esto que se revisa la puntuación MAPE del tercer modelo:

El resultado es el siguiente:

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.6925767

De acuerdo a los resultados del tercer modelo, estos presentan una pequeña mejora en el conjunto de entrenamiento y el conjunto de prueba, con un 2,81% y 5,07%, respectivamente.

Selección de modelos

Lo siguiente es seleccionar un modelo. De acuerdo a la puntuaciòn que obtuvo cada modelo, se conoce que el segundo y el tercer modelo poseen mejor rendimiento que el primero, sin embargo, el rendimiento entre estos no dista de mucho, ya que su puntuaciòn fue similar, con la diferencia de que el tercer modelo posee una pequeña ventaja de 0.2%. De acuerdo a estos resultados es pertinente cuestionarse si vale la pena utilizar la variable de retardo (es decir, la pérdida de 365 observaciones y el coste adicional de un grado de libertad para el modelo) para una mejora de menos del 0.2% Para esta pregunta no hay una respuesta directa, además de que dependerá de la previsión. Se tendrán en cuenta las siguientes pruebas:

El resultado es el siguiente:

summary(md_tslm3)$coefficients %>% tail(1)
##           Estimate Std. Error   t value     Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05

En este caso, el valor p de la variable lag365 quiere decir que esta es estadísticamente significativa a un nivel de 0.001.

El resultado es el siguiente:

anova(md_tslm3)

La prueba ANOVA también indica que la variable lag365 es estadísticamente significativa.

Para seleccionar las características se pueden aplicar métodos más fuertes, por ejemplo, la regresión por pasos, la cresta o el lazo.

Sin embargo, por simplicidad, para pronosticar la serie se utilizará el criterio de precisión y se seleccionará el tercer modelo. Lo siguiente será volver a entrenar el modelo en todas las series y pronosticar los siguientes 365 días:

final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365,
data = UKdaily)

Análisis de los residuos

Justo antes de finalizar la previsión, vamos a realizar un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals:

Obtenemos el siguiente resultado:

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-16

De acuerdo a lo que se observa en el gráfico de residuos, estos no son tumulto blanco, dado que existe cierta autocorrelación entre las series de residuos y sus rezagos. Lo anterior se interpreta como una indicación de que el modelo no capturó todos los patrones o la información que existe en las series. Una de las maneras en que se puede solucionar es identificando variables que puedan explicar la variación de los residuos, sin embargo, uno de los principales problemas es la dificultad que representa identificar las variables externas que pueden explicar el comportamiento de los residuos, y que a la vez, sean factibles de pronosticar.

En el caso en los que los patrones sobran en los residuos del modelo, entonces los residuos del modelo se deben tratar como datos de una serie temporal independiente, y además, modelizarse con otro modelo de previsión de series temporales. De igual manera se le puede dar un enfoque común, con una regresión con error ARIMA.

Finalización de la previsión

Por último, para finalizar el proceso y usar el modelo entrenado para predecir las próximas 365 observaciones. Como anteriormente ya se han utilizado variables externas como la función tslm, se tendrán que crear valores futuros. Se creará un objeto data.frame con los valores wday, month y lag365 para las siguientes 365 futuras observaciones. Este enfoque simplista involucra que primero se crean lo que corresponde a las observaciones previstas, y luego, extraer 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))

Lo siguiente es usar la variable date para crear las variables wday y month 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)

Se utiliza la variable forecast para crear la previsión:

UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)

Por último, se traza la previsión final con la función plot_forecast del paquete paquete TSstudio:

El resultado es el siguiente:

plot_forecast(UKgrid_fc,
title = "The UK National Demand for Electricity Forecast",
Ytitle = "MW",
Xtitle = "Year")

Parte 2

dataterremoto <- read.csv("database.csv")
class(dataterremoto)
## [1] "data.frame"
library(dplyr)
dataterremoto <- dataterremoto  %>%
  arrange(Date)
dataterremoto <- ts(dataterremoto$Magnitude, start = c(1965,1), end = c(2016,12), frequency = 12)
class(dataterremoto)
## [1] "ts"
library(TSstudio)
data(dataterremoto)
## Warning in data(dataterremoto): data set 'dataterremoto' not found
ts_plot(dataterremoto,
        title = "Terremotos entre 1965-2016 ocurridos",
        Ytitle = "Magnitud",
        Xtitle = "Fecha")

Esta es la representación gráfica de los terremotos ocurridos entre el año 1965 y el 2016. Como se puede observar, no hay una frecuencia directa entre el intervalo utilizado, por lo tanto no se puede asignar un número exacto de ocurrencia.

ts_info(dataterremoto)
##  The dataterremoto series is a ts object with 1 variable and 624 observations
##  Frequency: 12 
##  Start time: 1965 1 
##  End time: 2016 12
ts_decompose(dataterremoto)

En el grafico anterior, podemos observar que la serie presenta una tendencia bastante plana en los valores finales de 1990 e inicios del 2000. En este caso, igual que el anterior, los datos presentados no presentan una frecuencia que nos permita conocer la periocidad de los terremotos

dataterremoto_df <- ts_to_prophet(dataterremoto)
head(dataterremoto_df)
dataterremoto_df$trend <- 1:nrow(dataterremoto_df)
library(lubridate)
dataterremoto_df$seasonal <- month(dataterremoto_df$ds, label = T)
head(dataterremoto_df)
h <- 12 
train <- dataterremoto_df[1:(nrow(dataterremoto_df) - h), ]
test <- dataterremoto_df[(nrow(dataterremoto_df) - h + 1):nrow(dataterremoto_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.4523 -0.3027 -0.1563  0.1497  1.9512 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9541963  0.0357910 166.360   <2e-16 ***
## trend       -0.0001675  0.0001012  -1.656   0.0982 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4422 on 610 degrees of freedom
## Multiple R-squared:  0.004476,   Adjusted R-squared:  0.002844 
## F-statistic: 2.743 on 1 and 610 DF,  p-value: 0.09822
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 = "Fecha"),
           yaxis = list(title = "Magnitud"),
           legend = list(x = 0.05, y = 0.95))
  return(p)
}
plot_lm(data = dataterremoto_df, 
        train = train, 
        test = test,
        title = "Predicción del comportamiento de los terremotos. Magnitud vs. Tiempo")

En general, sabemos que este modelo es capaz de percibir el movimiento de la tendencia. Sin embargo, en esta tendencia lineal no se muestra una ruptura estructural como fue presentando en la primera parte, la línea mantiene su continuidad a lo largo de los años, disminuyendo en magnitud a medida que la fecha aumenta.

mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
                mean(abs(test$y - test$yhat) / test$y))
mape_trend
## [1] 0.05435221 0.04106426
md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)
## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4980 -0.3054 -0.1569  0.1657  1.9588 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.902843   0.017852 330.646   <2e-16 ***
## seasonal.L  -0.067736   0.061843  -1.095   0.2738    
## seasonal.Q  -0.138483   0.061843  -2.239   0.0255 *  
## seasonal.C   0.015539   0.061843   0.251   0.8017    
## seasonal^4   0.001288   0.061843   0.021   0.9834    
## seasonal^5  -0.039538   0.061843  -0.639   0.5228    
## seasonal^6   0.016671   0.061843   0.270   0.7876    
## seasonal^7   0.037590   0.061843   0.608   0.5435    
## seasonal^8   0.118299   0.061843   1.913   0.0562 .  
## seasonal^9   0.028090   0.061843   0.454   0.6498    
## seasonal^10  0.091641   0.061843   1.482   0.1389    
## seasonal^11  0.061945   0.061843   1.002   0.3169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4416 on 600 degrees of freedom
## Multiple R-squared:  0.02311,    Adjusted R-squared:  0.005199 
## F-statistic:  1.29 on 11 and 600 DF,  p-value: 0.2257
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)
plot_lm(data = dataterremoto_df, 
        train = train, 
        test = test,
        title = "Predicción del comportamiento de los terremotos. Magnitud vs. Tiempo")

En esta gráfica se captura la estructura del patrón estacional de la serie anteriormente presentada, como se puede observar, se ajusta de tal manera que sigue el patrón señalado por la magnitud. Sin embargo, la línea es recta, lo cual no se ajusto totalmente al comportamiento real.

mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
                   mean(abs(test$y - test$yhat) / test$y))
mape_seasonal
## [1] 0.05383805 0.04855883
md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)
## 
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5457 -0.2986 -0.1457  0.1629  1.9132 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9535528  0.0357037 166.749   <2e-16 ***
## seasonal.L  -0.0657571  0.0617677  -1.065   0.2875    
## seasonal.Q  -0.1384831  0.0617559  -2.242   0.0253 *  
## seasonal.C   0.0155388  0.0617559   0.252   0.8014    
## seasonal^4   0.0012884  0.0617559   0.021   0.9834    
## seasonal^5  -0.0395381  0.0617559  -0.640   0.5223    
## seasonal^6   0.0166715  0.0617559   0.270   0.7873    
## seasonal^7   0.0375903  0.0617559   0.609   0.5430    
## seasonal^8   0.1182989  0.0617559   1.916   0.0559 .  
## seasonal^9   0.0280899  0.0617559   0.455   0.6494    
## seasonal^10  0.0916410  0.0617559   1.484   0.1384    
## seasonal^11  0.0619448  0.0617559   1.003   0.3162    
## trend       -0.0001654  0.0001009  -1.639   0.1017    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.441 on 599 degrees of freedom
## Multiple R-squared:  0.02747,    Adjusted R-squared:  0.007988 
## F-statistic:  1.41 on 12 and 599 DF,  p-value: 0.1563
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = dataterremoto_df, 
        train = train, 
        test = test,
        title = "Predicción del comportamiento de los terremotos. Magnitud vs. Tiempo")

Teniendo en cuenta las gráficas anteriores, el siguiente paso es unir los dos componentes en un modelo y pronosticar los valores característicos de la serie. De esta manera, siguiendo la línea presentada en la primera gráfica, se pubo obtenar esta.

mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
              mean(abs(test$y - test$yhat) / test$y))
mape_md1
## [1] 0.05375073 0.04350817
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.6391 -0.2895 -0.1362  0.1768  1.8339 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.066e+00  5.335e-02 113.704  < 2e-16 ***
## seasonal.L  -6.576e-02  6.141e-02  -1.071  0.28470    
## seasonal.Q  -1.385e-01  6.140e-02  -2.257  0.02440 *  
## seasonal.C   1.554e-02  6.140e-02   0.253  0.80029    
## seasonal^4   1.288e-03  6.140e-02   0.021  0.98327    
## seasonal^5  -3.954e-02  6.140e-02  -0.644  0.51985    
## seasonal^6   1.667e-02  6.140e-02   0.272  0.78608    
## seasonal^7   3.759e-02  6.140e-02   0.612  0.54062    
## seasonal^8   1.183e-01  6.140e-02   1.927  0.05449 .  
## seasonal^9   2.809e-02  6.140e-02   0.458  0.64748    
## seasonal^10  9.164e-02  6.140e-02   1.493  0.13608    
## seasonal^11  6.194e-02  6.140e-02   1.009  0.31343    
## trend       -1.266e-03  4.019e-04  -3.149  0.00172 ** 
## I(trend^2)   1.795e-06  6.349e-07   2.827  0.00486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4385 on 598 degrees of freedom
## Multiple R-squared:  0.0403, Adjusted R-squared:  0.01943 
## F-statistic: 1.932 on 13 and 598 DF,  p-value: 0.02436
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = dataterremoto_df, 
        train = train, 
        test = test,
        title = "Predicción del comportamiento de los terremotos. Magnitud vs. Tiempo")

A partir de este momento, se ha añadido un ajuste polinomial a la gráfica, lo cual nos permité gráficar con mejor precisión el comportamiento de los terremotos, es decir, el modelo de regresión. Teniendo en cuenta que este ajuste reduce en el error en el conjunto de pruebas.

mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
              mean(abs(test$y - test$yhat) / test$y))
mape_md2
## [1] 0.05325876 0.05659700
dataterremoto_split <- ts_split(dataterremoto, sample.out = h)
train.ts <- dataterremoto_split$train
test.ts <- dataterremoto_split$test
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.6391 -0.2895 -0.1362  0.1768  1.8339 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.035e+00  7.917e-02  76.231  < 2e-16 ***
## season2     -1.550e-02  8.683e-02  -0.179  0.85836    
## season3      9.252e-02  8.683e-02   1.066  0.28707    
## season4      3.976e-02  8.683e-02   0.458  0.64723    
## season5      8.973e-02  8.683e-02   1.033  0.30184    
## season6      1.264e-01  8.683e-02   1.455  0.14610    
## season7      3.830e-02  8.683e-02   0.441  0.65930    
## season8      6.984e-02  8.683e-02   0.804  0.42157    
## season9     -5.746e-02  8.683e-02  -0.662  0.50843    
## season10     1.172e-01  8.684e-02   1.350  0.17760    
## season11    -6.696e-02  8.684e-02  -0.771  0.44098    
## season12    -5.896e-02  8.684e-02  -0.679  0.49739    
## trend       -1.266e-03  4.019e-04  -3.149  0.00172 ** 
## I(trend^2)   1.795e-06  6.349e-07   2.827  0.00486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4385 on 598 degrees of freedom
## Multiple R-squared:  0.0403, Adjusted R-squared:  0.01943 
## F-statistic: 1.932 on 13 and 598 DF,  p-value: 0.02436
r <- which(dataterremoto_df$ds == as.Date("2014-01-01"))
dataterremoto_df$s_break <- ifelse(year(dataterremoto_df$ds) >= 2010, 1, 0)
dataterremoto_df$s_break[r] <- 1
md3 <- tslm(dataterremoto ~ season + trend + I(trend^2) + s_break, data = dataterremoto_df)
summary(md3)
## 
## Call:
## tslm(formula = dataterremoto ~ season + trend + I(trend^2) + 
##     s_break, data = dataterremoto_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6397 -0.2919 -0.1320  0.1736  1.8357 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.040e+00  8.040e-02  75.128  < 2e-16 ***
## season2     -1.143e-02  8.563e-02  -0.134  0.89383    
## season3      9.059e-02  8.563e-02   1.058  0.29052    
## season4      5.030e-02  8.563e-02   0.587  0.55714    
## season5      8.578e-02  8.563e-02   1.002  0.31688    
## season6      1.236e-01  8.563e-02   1.443  0.14957    
## season7      4.288e-02  8.564e-02   0.501  0.61678    
## season8      6.603e-02  8.564e-02   0.771  0.44096    
## season9     -5.505e-02  8.564e-02  -0.643  0.52062    
## season10     1.162e-01  8.564e-02   1.357  0.17541    
## season11    -5.683e-02  8.564e-02  -0.664  0.50723    
## season12    -6.061e-02  8.565e-02  -0.708  0.47940    
## trend       -1.373e-03  4.687e-04  -2.930  0.00352 ** 
## I(trend^2)   2.062e-06  8.338e-07   2.473  0.01365 *  
## s_break     -9.285e-02  8.792e-02  -1.056  0.29135    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4366 on 609 degrees of freedom
## Multiple R-squared:  0.03863,    Adjusted R-squared:  0.01653 
## F-statistic: 1.748 on 14 and 609 DF,  p-value: 0.04302
ts_heatmap (dataterremoto ,  last  =  NULL ,  wday  =  TRUE ,  color  =  "Blues" , 
  title  = "Terremotos entre 1965-2016"  ,  padding  =  TRUE )

Finalmente, esta gráfica presenta los datos estaciones, clasificados según su mes teniendo en cuenta su año. A su vez, cada recuadro contiene la magnitud que se obtuvo en ese periodo específico.

library(dplyr)
dataterremoto_df<- dataterremoto_df%>%
  mutate(wday = wday(ds, label = TRUE),
         month = month(ds, label = TRUE),
         lag365 = dplyr::lag(y   , 365)) %>%
  filter(!is.na(lag365)) %>%
  arrange(ds)
str(dataterremoto_df)
## 'data.frame':    259 obs. of  8 variables:
##  $ ds      : Date, format: "1995-06-01" "1995-07-01" ...
##  $ y       : num  5.6 6 5.7 5.7 6 6 5.8 6 5.5 5.8 ...
##  $ trend   : int  366 367 368 369 370 371 372 373 374 375 ...
##  $ seasonal: Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 6 7 8 9 10 11 12 1 2 3 ...
##  $ s_break : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ wday    : Ord.factor w/ 7 levels "dom\\."<"lun\\."<..: 5 7 3 6 1 4 6 2 5 6 ...
##  $ month   : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 6 7 8 9 10 11 12 1 2 3 ...
##  $ lag365  : num  6.5 5.9 5.6 5.6 6 6.8 6 5.9 5.7 6.2 ...
start_date <- min(dataterremoto_df$ds)
start <- c(year(start_date), yday(start_date))
dataterremoto_ts <- ts(dataterremoto_df$y,
 start = start,
 frequency = 365)
acf(dataterremoto_ts, lag.max = 365 * 51)

h <- 12
dataterremoto_partitions <- ts_split(dataterremoto_ts, sample.out = h)
train_ts <- dataterremoto_partitions$train
test_ts <- dataterremoto_partitions$test
train_df <- dataterremoto_df[1:(nrow(dataterremoto_df) - h), ]
test_df <- dataterremoto_df[(nrow(dataterremoto_df) - h + 1):nrow(dataterremoto_df), ]
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365,
data = train_df)
summary(md_tslm3)$coefficients %>% tail(1)
##           Estimate Std. Error t value Pr(>|t|)
## season365      1.3        NaN     NaN      NaN
anova(md_tslm3)
## Warning in anova.lm(md_tslm3): ANOVA F-tests on an essentially perfect fit are
## unreliable
final <- tslm(dataterremoto_ts ~ season + trend + wday + month + lag365,
 data = dataterremoto_df)
checkresiduals(dataterremoto)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Conclusión

En el informe realizado, cuyo énfasis fue la aplicación de previsión del modelo de regresión lineal en el software Rstudio, se emplearon una gran variedad de métodos con el fin de generar conocimiento para cuantificar y reconocer la relación que existe entre las variables dependientes, únicas e independientes múltiples ya que fueron de gran ayuda para ponerla en practica al momento de realizar los diferentes tipos de gráficas, los datos con los que se implementaron las graficas en la Parte 2 se basan en una toma de datos sobre terremotos ocurridos durante un intervalo de años, que va desde 1965 a 2019, todo estos datos (tanto los de la Parte 1 y la Parte 2) fueron aplicados de manera clara y correcta en el Software Rstudio y serán puestos en practica para generar los datos que se requieren en cada incisos del informe. Por último, tras finalizar el desarrrollo de ambas partes del informe se reconoce la importancia de implementar un modelo de regresión lineal, así como las entradas de regresión lineal que tienen como objetivo representar la relación entre dos o más variables, que a su vez permitieron realizar predicciones de los valores que estas tomaran.