Parte 1

El modelo de regresión lineal es uno de los métodos más comunes para identificar y cuantificar la relación entre una variable dependiente y una sola (regresión lineal univariante) o múltiples (regresión lineal univariante) o múltiples (regresión lineal multivariante) variables independientes. Este modelo tiene una amplia gama de aplicaciones, desde la inferencia causal hasta el análisis predictivo y, en en particular, la previsión de series temporales.

Este capítulo se centra en los métodos y enfoques de previsión de datos de series temporales con regresión lineal. Esto incluye métodos para descomponer y pronosticar los componentes de la serie de la serie (por ejemplo, la tendencia y los patrones estacionales), el tratamiento de eventos especiales (como como los valores atípicos y los días festivos), y el uso de variables externas como regresores

Este capítulo abarca los siguientes temas:

$\cdot$ Enfoques de previsión con modelos de regresión linea

$\cdot$ Extracción y estimación de los componentes de la serie

$\cdot$ Tratamiento de las rupturas estructurales, los valores atípicos y los eventos especiales

$\cdot$ Previsión de series con multiestacionalidad

Requisitos técnicos

En este capítulo se utilizarán los siguientes paquetes:

$\cdot$ Tstudio: Versión 0.1.4 y superior

$\cdot$ plotly: Versión 4.8 y superior

$\cdot$ dplyr: Versión 0.8.1 y superior

$\cdot$ lubridate: Versión 1.7.4 y superior

$\cdot$ forecast: Versión 8.5 y superior

Regresión lineal

El uso principal del modelo de regresión lineal es cuantificar la relación entre la variable dependiente Y (también conocida como variable de respuesta) y la/s variable/s independiente/s independientes X (también conocidas como variables predictoras, impulsoras o regresoras) de forma lineal. En otras palabras, el modelo expresa la variable dependiente como una combinación lineal de las variables independientes. Una relación lineal entre las variables dependientes e independientes puede generalizarse mediante las siguientes ecuaciones:

$\cdot$ En el caso de una sola variable independiente, la ecuación es la siguiente :

\[ \begin{aligned}Y_i = \beta_0 + \beta_1 * X_{1,i} + \varepsilon_i\end{aligned} \] $\cdot$ Para $n$ variables independientes, la ecuación tiene el siguiente aspecto:

\[ \begin{aligned}Y_i = \beta_0 + \beta_1X_{1,i} + \beta_2X_{2,i} + ... + \beta_nX_{n,i} + \varepsilon_i\end{aligned} \] Las variables del modelo para estas ecuaciones son las siguientes:

$\cdot$ $i$ representa el índice de observaciones $i = 1,..., N$

$\cdot$ $Y_i$ representa la observación i de la variable dependiente $i = 1,..., N$

$\cdot$ $X_{j,i}$ representa el valor $i$ de la variable independiente $j$, donde $j = 1,..., n$

$\cdot$ $\beta_0$ representa el valor del término constante (o intercepción)

$\cdot$ $\beta_j$ representa los parámetros correspondientes (o coeficientes) de las $j$ variables independientes

$\cdot$ $\varepsilon_i$ define el término de error, que no es más que toda la información que no fue capturada por las variables independientes para la observación $i$

El término lineal, en el contexto de la regresión, se refiere a los coeficientes del modelo, que deben seguir una estructura lineal (ya que esto nos permite construir una combinación lineal a partir de las variables independientes). Por otra parte, las variables independientes pueden seguir una formación tanto lineal como no lineal.

Suponiendo que las ecuaciones anteriores representen la verdadera naturaleza de la relación lineal entre las variables dependientes e independientes, entonces el modelo de regresión lineal proporciona una estimación de esos coeficientes (es decir, $\beta_0,\beta_1,...,\beta_n$), que se pueden formalizar mediante las siguientes ecuaciones:

$\cdot$ Para el modelo de regresión lineal univariante, la ecuación es la siguiente:

\[ \begin{aligned}\hat{Y}_i = \hat{\beta_0} + \hat{\beta_1} X_{1,i} \end{aligned} \] $\cdot$ Para el modelo de regresión lineal multivariante, la ecuación es la siguiente:

\[ \begin{aligned}\hat{\varepsilon_i} = \hat{Y_i} + \hat{\beta_1}X_{1,i} + \hat{\beta_2}X_{2,i} + ... + \hat{\beta_n}X_{n,i} \end{aligned} \] Las variables del modelo para estas ecuaciones son las siguientes:

$\cdot$ $i$ representa el índice de observaciones $i = 1,..., N$

$\cdot$ $\hat{Y_i}$ representa la estimación de la observación de la variable dependiente $i$

$\cdot$ $X_{j,i}$ representa el valor $i$ de la variable independiente $j$, donde $j = 1,..., n$

$\cdot$ $\hat{\beta_0}$ representa la estimación del término constante (o intercepción)

$\cdot$ $\hat{\beta_1},...,\hat{\beta_n}$ son la estimación de los correspondientes parámetros (o coeficientes) de las $n$ variables independientes

La estimación de los coeficientes del modelo se basa en los dos pasos siguientes:

$\cdot$ Definir una función de coste (también conocida como función de pérdida)-estableciendo alguna métrica de error para minimizar

$\cdot$ Aplicar la optimización matemática para minimizar la función de coste

El enfoque de estimación más común es la aplicación de los mínimos cuadrados ordinarios (OLS) como técnica de optimización para minimizar la suma de cuadrados de los residuos \[ \sum_{i=1}^{N}\hat{\varepsilon_i^2} \]

Elevar al cuadrado los residuos tiene dos efectos:

$\cdot$ Evita que los valores positivos y negativos se anulen entre sí al sumarlos los valores positivos y negativos se cancelen entre sí al sumarlos.

$\cdot$ El efecto cuadrado proporciona una penalización exponencial para los residuos con mayor distancia, ya que su coste es mayor

Existen múltiples técnicas de estimación, además de las OLS, como la de máxima probabilidad, el método de los momentos y el bayesiano. Aunque esos métodos no entran en el ámbito del libro, es muy recomendable leer y conocer los enfoques alternativos

Estimación de los coeficientes con el método OLS

El OLS es un método de optimización sencillo que se basa en el álgebra lineal y el cálculo básicos, o cálculo matricial. (Esta sección es de conocimiento general; si no está familiarizado con el cálculo matricial puede saltarse esta sección). El objetivo del OLS es identificar los coeficientes que minimicen la suma de cuadrados de los residuos. Supongamos que el residuo de la observación $i$ se define de la siguiente manera:

\[ \begin{aligned}\hat{\varepsilon_i} =Y_1 -\hat{Y_1} \end{aligned} \] Podemos entonces fijar la función de costes con la siguiente expresión: \[\begin{aligned} \sum_{i=1}^{N}\hat{\varepsilon_i^2}=(Y_1 -\hat{Y_1})^2+(Y_2-\hat{Y_2})^2+...+(Y_n -\hat{Y_n})^2\end{aligned} \] Antes de aplicar el método OLS para minimizar la suma de cuadrados de los residuos, por simplicidad, transformaremos el representante de la función de costes en una matriz de formación:

$$ \[\begin{equation*} Y= \begin{bmatrix} Y_1 \\ Y_2 \\ \cdot \\ \cdot \\ \cdot \\ Y_N\\ \end{bmatrix} X=\begin{bmatrix}1 & X_{1,1} &\cdot& \cdot& \cdot& X_{1,n}\\ \cdot \\ \cdot \\ \cdot \\ 1& X_{N,1} &\cdot& \cdot& \cdot& X_{N,n}\\ \end{bmatrix} \beta= \begin{bmatrix} \beta_0\\ \beta_1 \\ \cdot \\ \cdot \\ \cdot \\ \beta_n\ \end{bmatrix} \varepsilon=Y-X\beta= \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \cdot \\ \cdot \\ \cdot \\ \varepsilon_N\\ \end{bmatrix} \end{equation*}\]$$

Estos conjuntos de matrices representan lo siguiente:

$\cdot$ Vector Y (o $N \times 1$ matriz), que representa una variable dependiente con $N$ observaciones

$\cdot$ X, una matriz de dimensiones $N \times n+1$, que representa las correspondientes $n$ variables independientes y un escalar de 1’s para el componente de intercepción ($\beta_0$)

$\cdot$ $\beta$,una matriz de dimensiones $n+1 \times 1$, que representa los coeficientes del modelo

$\cdot$ $\varepsilon$,una matriz de dimensiones $N \times 1$, que representa el error correspondiente (o la diferencia entre ellos) del valor actual de $Y$ y estimar $\beta X$

El término residual $ $ no debería ser confundida por error con el término ${\epsilon_i}$. La primera representa las diferencias entre las series $Y$ y el estimado de $\hat{Y}$; por otro lado, el segundo término (término de error) representa las diferencias entre las series y el valor esperado.

Se establece la función de costo usando la forma matricial la cual se definió anteriormente:

\[ \sum{\epsilon^2}=\sum{\epsilon^\tau\epsilon} \]

Se prosigue a expandir esta expresión a partir de la fórmula de $\epsilon$:

\[ {\epsilon^\tau\epsilon} =(Y-X\beta)^T(Y-X\beta) \]

Asimismo, se multiplican los dos componentes $(\epsilon^\tau$ y $\epsilon)$ y se abren los paréntesis:

\[ {\epsilon^\tau\epsilon} =Y^TY-2Y^TX\beta+\beta^TX^TX\beta \]

Dado que el objetivo es encontrar el $\beta$ que minimiza esta ecuación, se diferencia la ecuación con respecto a $\beta$ y luego se iguala a 0:

\[ \left(\frac{\partial{\epsilon^\tau\epsilon}}{\partial\beta}\right) = \left(\frac{\partial(Y^TY-2Y^TX\beta+\beta^TX^TX\beta)}{\partial\beta}\right) = 0 \]

Al resolver la ecuación, se obtiene lo siguiente,

\[ X^TX\beta=X^TY \] Por otro lado, al manipular esta ecuación permite extraer la matriz, la estimación de la matriz de coeficientes:

\[ \hat{\beta} = (X^TX)^{-1}X^TY \] Es importante notar que se cambia la notación de $\beta$ a $\hat{\beta}$ para representar el verdadero valor estimado de $\beta$.

Las principales propiedades de los coeficientes de estimación OLS son:

$\cdot$ La característica principal del método de estimación de los coeficientes de OLS es la falta de sesgo de la estimación de los coeficientes con respecto a los valores reales. En otras palabras, para cualquier, $\hat{\beta_i}, E(\hat{\beta_i})= {\beta_i}$

$\cdot$ La línea de regresión muestral siempre pasará por la media de $X$ e $Y$

$\cdot$ El valor de $\hat{Y}$ la estimación OLS de la variable dependiente $Y$, es igual a $\bar{Y}$, la media la variable dependiente

$\cdot$ La media del vector $\hat{\epsilon}$ de residuos es igual a cero, o \[\frac{\sum_{i=1}^{N}\hat{\epsilon_i}}{N}\]

Supuestos del OLS

El modelo de OLS tiene los siguientes supuestos:

$\cdot$ El modelo de los coeficientes deben seguir una estructura lineal (por ejemplo, $Y=\hat{\beta_0} + \hat{\beta_1}e^x$ es un modelo lineal pero $Y=\hat{\beta_0} + X_1^{\hat{\beta_1}}$)

$\cdot$ No existe una colinealidad perfecta entre variables independientes $X_1,X_2,..,X_n$. En otras palabras, ninguna de las variables independientes es una combinación lineal de ninguna de las otras variables independientes.

$\cdot$ Todas las variables independientes deben tener una varianza distinta de cero (o no constante).

$\cdot$ El término de error$\epsilon$, condicionado a la matriz de variables $X$ independientes X, es una variable independiente e idénticamente distribuida (i.i.d) con media 0 y varianza constante $\sigma^2$

$\cdot$ Tanto las variables dependientes como las independientes se extraen de la población en una muestra aleatoria. Este supuesto no se cumple cuando se hacen regresiones de datos de series de tiempo, ya que típicamente las observaciones tienen cierto grado de correlación. Por lo tanto, esta suposición se relaja al hacer una regresión de datos de series de tiempo.

Pronóstico con regresión lineal

El modelo de regresión lineal, no fué diseñado para manejar ni pronosticar datos de series de tiempo. Caso contrario sucede con los modelos tradicionales de series de tiempo como ARIMA o Holt-Winters. En cambio, es un modelo genérico con muchas aplicaciones, las cuales van desde la inferencia causal hasta el análisis predictivo.

Por lo tanto, la previsión con un modelo de regresión lineal se basa principalmente en:

$1.$ Identificar la estructura, características clave, patrones, valores atípicos de la serie.

$2.$ Transformar esas características en variables de entrada y hacer una regresión con la serie para crear un modelo de pronóstico.

Las características principales de un modelo de pronóstico de regresión lineal son los componentes de tendencia y estacionales.

Pronosticar la tendencia y los componentes estacionales

En el Capítulo 5, Descomposición de datos de series de tiempo, se presentaron los componentes estructurales de la serie, los cuales son: la tendencia, el ciclo, los componentes estacionales e irregulares y los métodos para descomponerlos con la función de descomposición.

La descomposición de una serie se puede definir mediante:

$\cdot$ ${Y_i}= T_t + S_t + C_t + I_ t$, sí la serie tiene una estructura aditiva.

$\cdot$ ${Y_i}= T_t \times S_t \times C_t \times I_ t$, sí la serie tiene una estructura multiplicativa.

Se tiene que,

Tendencia: representa el crecimiento de la serie a lo largo del tiempo después de ajustar y eliminar los efectos estacionales.

Estacional: representa un patrón cíclico recurrente el cual se deriva directamente de las unidades de frecuencia de la serie (por ejemplo, el mes del año para una serie con una frecuencia mensual).

Ciclo: es un patrón cíclico que no está relacionado con la unidad de frecuencia de la serie.

Irregular: representa cualquier otro patrón que no sea capturado por los componentes demás componente (tendencia, estacional y ciclo).

Para mayor simplicidad, se puede descartar el componente de ciclo, ya que normalmente se fusionó con el componente de tendencia (o se ignoró el componente de ciclo). Por lo cual, se puede actualizar y reemplazar las ecuaciones anteriores con lo siguiente:

$\cdot$ ${Y_i}= T_t + S_t + I_t$, sí la serie tiene una estructura aditiva.

$\cdot$ ${Y_i}= T_t \times S_t \times I_t$, sí la serie tiene una estructura multiplicativa.

Asimismo, se pueden transformar esas ecuaciones para un modelo de regresión lineal, en donde se cambian la notación de las ecuaciones,

\[ Y_t ={\beta_0} + {\beta_1T_t} + {\beta_2S_t} + {\epsilon_t} \] donde,

$\cdot$ $Y$ representa una serie de tiempo con n observaciones.

$\cdot$ $T$ es una variable independiente con n observaciones que representa el componente de tendencia de la serie.

$\cdot$ $S$ es una variable independiente con n observaciones y representa el componente estacional de la serie

$\cdot$ $\epsilon_t$ es el término de error de la regresión.Representa el componente irregular o cualquier patrón que no sea tomado por la tendencia de la serie y el componente estacional.

$\cdot$ $\beta_0,\beta_1$ y $\beta_2$ son aquellos que representan la intersección del modelo y los coeficientes de los componentes de tendencia y estacional, respectivamente.

Nota: Se cambia la notación de observaciones de i a t, ya que esta representa la dimensión de tiempo de la serie.

Para transformar una serie con una estructura multiplicativa en una formación de regresión lineal se tuvo necesario contar con una transformación de la serie en una estructura aditiva. Esto se puede lograr al aplicar la transformación logarítmica para ambos lados de las ecuaciones:

\[ Log(Y_t) ={Log(T_t) + Log(S_t) + Log(I_t)} \]

Posterior a que la serie se transformó en una estructura aditiva, para transformar a una formación de regresión lineal se continúa el mismo proceso que se estudió anteriormente:

\[ Log(Y_t) ={\beta_0}+{\beta_1Log(T_t) + \beta_2Log(S_t) + Log(I_t)} + {\epsilon_t} \]

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

En los siguientes ejemplos, se muestra el proceso de creación de nuevas funciones a partir de series existentes utilizando la serie USgas.

Se obtiene la serie desde el paquete TSstudio nuevamente y trazémosla con la función ts_plot:

library(TSstudio)

data("USgas")

ts_plot(USgas,
        title = "US Monthly Natural Gas consumption",
        Ytitle = "Billion Cubic Feet",
        Xtitle = "Year")

Además, se revisan las principales características de la serie haciendo uso de 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 10

Como se puede observar en el gráfico de la serie,USgas es una serie mensual con un fuerte componente estacional mensual y una línea de tendencia bastante estable. Asimismo, se puede estudiar más la estructura de los componentes de la serie con la función ts_decompose:

ts_decompose(USgas)

Al observar el gráfico se puede denotar que la tendencia de la serie es bastante plana entre 2000 y 2010, y en el futuro, tiene un crecimiento lineal. Por lo cual, la tendencia general entre 2000 y 2018 no es necesariamente lineal.

Antes de usar la función lm, la función de regresión lineal R incorporada del paquete stats, se tiene que transformar la serie de un objeto ts a un objeto data.frame. Lo anterior permite inferir que se hará uso de la función ts_to_prophet del paquete TSstudio:

USgas_df <- ts_to_prophet(USgas)

La función transforma el objeto ts dentro de dos columnas de data.frame, en la cual, una de las columnas representan el tiempo y la otra columna los componentes de la serie:

"Se obtiene la siguiente salida: "

## [1] "Se obtiene la siguiente salida: "

head(USgas_df)

##           ds      y
## 1 2000-01-01 2510.5
## 2 2000-02-01 2330.7
## 3 2000-03-01 2050.6
## 4 2000-04-01 1783.3
## 5 2000-05-01 1632.9
## 6 2000-06-01 1513.1

Después de transformar la serie en un objeto de data.frame, se puede empezar a crear las características de entrada de la regresión. La primera característica que se crea se le denomina Tendencia de la Serie. Un enfoque básico para construir la variable de tendencia es indexar de forma cronológica, las observaciones de la serie:

USgas_df$trend <- 1:nrow(USgas_df)

La regresión de la serie con su índice proporciona una estimación del crecimiento marginal de un mes a otro, ya que el índice se encuentra en un orden cronológico con un incremento constante.

La segunda característica de la regresión es el componente estacional. Debido a que se quiere medir la contribución de cada unidad de frecuencia a la oscilación de serie, se hace uso de una variable categórica para cada unidad de frecuencia. En el caso de la serie USgas, las unidades de frecuencia representan los meses del año, por lo cual se crea una variable categórica con 12 categorías, en donde cada categoría es un mes del año. Por otro lado, se usa la función de month del paquete Lubridate para extraer los meses 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 <- month(USgas_df$ds, label = T)

Asimismo, se hace uso de la función factorial para convertir la salida de la función de mes en una variable categórica sin orden. Se procede a revisar los datos después de agregar las funciones mencionadas anteriormente:

head(USgas_df)

##           ds      y trend seasonal
## 1 2000-01-01 2510.5     1      Jan
## 2 2000-02-01 2330.7     2      Feb
## 3 2000-03-01 2050.6     3      Mar
## 4 2000-04-01 1783.3     4      Apr
## 5 2000-05-01 1632.9     5      May
## 6 2000-06-01 1513.1     6      Jun

Por otro lado, antes de empezar a realizar la regresión de la serie con las características en mención, se divide la serie de tiempo en una partición de entrenamiento y prueba. Se configura los ultimos 12 meses de la serie como partición de prueba:

h <- 12 

train <- USgas_df[1:(nrow(USgas_df) - h), ]

test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]

Después de crear los marcos de datos de entrenamiento y prueba, se puede revisar como el modelo de regresión toma cada uno de los componentes juntos y separados.

Modelando la tendencia de la serie y estacional componentes

En primer lugar, se modela la tendencia de la serie, en la cual, se hace una regresión de la serie con la variable de tendencia, en la partición de entrenamiento:

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

Se hace uso de la función summary para revisar los detalles del modelo:

summary(md_trend)

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

Como se puede observar en el resultado anterior, el coeficiente de la variable de tendencia tiene un nivel de significancia de 0.001. Sin embargo, el R-Cuadrado ajustado de la regresión es bastante bajo, lo cual se relaciona con el patrón estacional que se observa en los gráficos anteriores.

Asimismo, se puede denotar en el resultado de la regresión anterior que la cuarta columna representa el valor-p de cada uno de los coeficientes del modelo. El valor-p brinda la probabilidad de que se rechace la hipótesis nula dada que es realmente verdadera, es decir, el error de tipo I. Por lo tanto, el el valor-p menor que $\alpha$, el valor umbral, se rechaza la hipótesis nula $H_o$ con un nivel de significancia de $\alpha$, donde los valores típicos de $\alpha$ son 0.1,0.05,0.01,etc.

Se prosigue a utilizar la función de forecast del paquete de estadísticos, el cual predice los valores de datos de entrada basados en un modelo determinado.

De igual manera, se usa md_trend para predecir los valores ajustados y previstos del modelo de tendencia anterior:

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

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

Se crea una función de utilidad que traza la serie y la salida del modelo, usando el paquete ploty:

library(plotly)

## 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 = ""),
           yaxis = list(title = "Billion Cubic Feet"),
           legend = list(x = 0.05, y = 0.95))
  return(p)
}

La función tiene los siguientes argumentos:

$\cdot$ data: los datos de entrada, un objeto data.frame que sigue la misma estructura que la de USgas (incluyendo la variable yhat ).

$\cdot$ train: Conjunto de entrenamiento correspondiente que se usó para entrenar el modelo.

$\cdot$ test: Conjunto de pruebas correspondiente que se utilizó para evaluar el modelo pronóstico.

$\cdot$ tittle: título de la trama, se utiliza por defecto NULL.

Se establece las entradas de la función plot_lm con el modelo de salida:

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

En general, el modelo fué capaz de tomar el movimiento general de tendencia, sin embargo, una tendencia lineal no puede capturar la ruptura estructural de la tendencia que ocurrió en el año 2010. Para el análisis de comparación, se desea medir la tasa de error del modelo de entrenamiento y conjuntos de prueba: La norma se mide por medio de:

$MAPE=(\frac{1}{M})\sum|\frac{A_i-F_i}{A_i}|$

mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
                mean(abs(test$y - test$yhat) / test$y))

"Se obtiene la siguiente salida: "

## [1] "Se obtiene la siguiente salida: "

mape_trend

## [1] 0.1644088 0.1299951

El proceso de modelar y pronosticar el componente estacional sigue el mismo proceso de tendencia, al realizar una regresión de la serie con la variable estacional creada anteriormente:

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

summary(md_seasonal)

## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -608.98 -162.34  -50.77  148.40  566.89 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2030.88      14.43 140.747  < 2e-16 ***
## seasonal.L   -480.00      50.24  -9.554  < 2e-16 ***
## seasonal.Q   1103.83      50.17  22.000  < 2e-16 ***
## seasonal.C     72.42      50.05   1.447 0.149389    
## seasonal^4    174.60      50.07   3.487 0.000592 ***
## seasonal^5    288.01      50.13   5.745 3.13e-08 ***
## seasonal^6    -44.67      50.09  -0.892 0.373460    
## seasonal^7   -187.91      49.96  -3.762 0.000218 ***
## seasonal^8     84.95      49.84   1.705 0.089706 .  
## seasonal^9     46.16      49.78   0.927 0.354828    
## seasonal^10    77.55      49.76   1.559 0.120587    
## seasonal^11   -11.09      49.75  -0.223 0.823856    
## ---
## 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

Debido a que se realizó una regresión de la variable dependiente con la categórica, el modelo de regresión crea coeficientes de 11 a 12 categorías, que son incluidas con los valores de la pendiente. De igual manera, todos los coeficientes del modelo son estadísticamente significativos. Por último, se puede observar que la función R-cuadrado, ajustado al modelo estacional es más alto que el modelo de tendencia (0.78 vs 0.1).

Se prosigue a actualizar 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)

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

Como se observa en el gráfico anterior, el modelo logra capturar adecuadamente la estructura del patrón estacional de la serie. Sin embargo, se puede denotar que hace falta la tendencia de la serie. Por lo tanto, antes de agregar la tendencia y los componentes estacionales, se estudia 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.21502100

Se puede decir que la alta tasa de error en el conjunto de pruebas se relaciona con el componente de tendencia que no se incluye en el modelo anterior. Por otro lado, el siguiente paso a realizar es unir los dos componentes en un modelo y posterior a esto, pronosticar los valores característicos de la serie:

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

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) 1733.7153    17.0794 101.509  < 2e-16 ***
## seasonal.L  -498.1709    29.6354 -16.810  < 2e-16 ***
## seasonal.Q  1115.2951    29.5872  37.695  < 2e-16 ***
## seasonal.C    78.9835    29.5103   2.676  0.00802 ** 
## seasonal^4   175.6505    29.5196   5.950 1.09e-08 ***
## seasonal^5   285.0192    29.5568   9.643  < 2e-16 ***
## seasonal^6   -49.3611    29.5319  -1.671  0.09610 .  
## seasonal^7  -192.3050    29.4540  -6.529 4.77e-10 ***
## seasonal^8    81.8787    29.3835   2.787  0.00581 ** 
## seasonal^9    44.4849    29.3480   1.516  0.13106    
## seasonal^10   76.8636    29.3372   2.620  0.00943 ** 
## seasonal^11  -11.2755    29.3353  -0.384  0.70109    
## 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 con los componentes de tendencia y estacional brinda un aumento adicional al R-cuadrado ajustado del modelo de 0,78 a 0,91. En el gráfico de la salida del modelo se denota lo anterior:

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


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

Se mide 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.09143290

La regresión de la serie con los componentes de tendencia y estacionales proveé un aumento significativo en la calidad de ajuste y en la precisión del modelo. A pesar de ello, al analizar la gráfica del ajuste y la previsión del modelo, se observa que la tendencia conserva una linealidad, así como también, falta la ruptura estructural de la tendencia de la serie.

Debido a lo anterior, para capturar una tendencia no lineal es necesario agregar un componente polinomial a la tendencia de la serie para obtener la curvatura de la tendencia a largo tiempo. Razón por la cual, se hace uso del argumento I para aplicar operaciones matemáticas en cualquiera de los dos objetos de entrada. Por lo tanto, se usa este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:

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

summary(md2)

## 
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.882e+03  2.199e+01  85.568  < 2e-16 ***
## seasonal.L  -4.917e+02  2.530e+01 -19.438  < 2e-16 ***
## seasonal.Q   1.121e+03  2.525e+01  44.381  < 2e-16 ***
## seasonal.C   8.247e+01  2.518e+01   3.275  0.00123 ** 
## seasonal^4   1.763e+02  2.519e+01   6.999 3.35e-11 ***
## seasonal^5   2.835e+02  2.522e+01  11.243  < 2e-16 ***
## seasonal^6  -5.175e+01  2.520e+01  -2.054  0.04123 *  
## seasonal^7  -1.946e+02  2.513e+01  -7.741 3.97e-13 ***
## seasonal^8   8.030e+01  2.507e+01   3.203  0.00157 ** 
## seasonal^9   4.362e+01  2.504e+01   1.742  0.08293 .  
## seasonal^10  7.651e+01  2.503e+01   3.057  0.00253 ** 
## seasonal^11 -1.137e+01  2.503e+01  -0.454  0.65005    
## 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

Al adicionar el polinomio de segundo grado al modelo de regresión no se llevo a cabo una mejora significativa de la bondad de ajuste del modelo. Por otro lado, puede observarse en la siguiente gráfica, este simple cambio en la estructura del modelo en el cual se obtiene 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 Seasonal Component of the Series")

En el siguiente modelo de puntuación MAPE, la precisión del modelo está significativamente mejorado al agregar la línea de tendencia polinómica al modelo de regresión, tal como el error en el conjunto de pruebas, el cual cayó de 9,2% a 4,5%, se obtiene lo siguiente:

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, se ha logrado ver el proceso manual de transformar un objeto ts en una regresión lineal de formato de modelo de predicción. La función tslm del paquete de previsión proporciona una función integrada con el objetivo de transformar un objeto ts en un modelo de predicción de regresión lineal. Haciendo uso de la función tslm, es posible configurar el componente de regresión junto con otras características.

Ahora se repetirá el ejemplo anterior y se pronosticarán las últimas 12 observaciones de la serie USgas con la función tslm utilizando la tendencia, el cuadrado de la tendencia, y los componentes estacionales. En primer lugar, se divide la serie en particiones de entrenamiento y pruebas usando la función ts_split:

USgas_split <- ts_split(USgas, sample.out = h)

train.ts <- USgas_split$train

test.ts <- USgas_split$test

A continuación, se aplicará la misma fórmula que se usó para crear el anterior modelo de predicción md2 usando la función tslm:

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

md3 <- tslm(train.ts ~ season + trend + I(trend^2))

Ahora se va a revisar el md3, el resultado de la función tslm, y compararla con el resultado de md2, obteniendo el 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

Como se puede observar en el resultado anterior, ambos modelos (md2 y md3) son exactamente iguales. Existen diversas ventajas con el uso de la función tslm, en lugar de configurar manualmente 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, además, ingeniería de características.
El objeto de salida soporta toda la funcionalidad del modelo de predicción (como las funciones de exactitud y checkresiduals) y paquetes de TSstudio (como las funciones test_forecast y plot_forecast)

Modelado de eventos individuales y eventos no estacionales

En algunas ocasiones, los datos de series temporales pueden o no contener patrones inusuales a través del tiempo. Los siguientes son ejemplos de tales eventos:

Valores atípicos: Evento(s) que están fuera de los patrones normales de la serie.

Ruptura estructural: es algún evento significativo que genera un cambio en los patrones históricos de la serie. Un ejemplo bastante común es un cambio en el crecimiento de la serie.
Eventos recurrentes no estacionales: sucede cuando un evento se repite de ciclo a ciclo, pero el momento en el que ocurren dichos cambios cambia de ciclo a ciclo. Un ejemplo común de este evento son las vacaciones de Pascua, las cuales ocurren de Marzo a Abril todos los años.

Al no aparecer en el modelo de regresión, este tipo de eventos sesgarán la estimación de los coeficientes, debido a que el modelo ponderará estos eventos junto con los eventos regulares con comportamiento regular de la serie. El uso del hot enconding, la codificación binaria, o variables de bandera podría ayudar al modelo, tanto para ignorar este tipo de eventos o ajustar los coeficientes del modelo.

Por ejemplo, es posible observar en la trama de descomposición de la serie USgas que fue expuesta anteriormente, que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Mientras que el crecimiento antes del 2010 fue relativamente nula, la pendiente de la tendencia cambió después, demostrando un crecimiento positivo. En este caso, es posible hacer uso de una variable binaria que sea igual a cero para las observaciones 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 con las variables correspondientes. El siguiente ejemplo expone el proceso de creación de una variable binaria externa (0) antes del año 2010 y 1 después, utilizando 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 se usará la nueva característica para remodelar la serie USgas:

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

Se usará la función de resumen para revisar la salida del modelo:

summary(md3)

## 
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break, 
##     data = USgas_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -469.25  -50.68   -2.66   63.63  275.89 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.661e+03  3.200e+01  83.164  < 2e-16 ***
## season2     -3.054e+02  3.448e+01  -8.858 2.61e-16 ***
## season3     -4.849e+02  3.448e+01 -14.062  < 2e-16 ***
## season4     -9.272e+02  3.449e+01 -26.885  < 2e-16 ***
## season5     -1.108e+03  3.449e+01 -32.114  < 2e-16 ***
## season6     -1.127e+03  3.450e+01 -32.660  < 2e-16 ***
## season7     -9.568e+02  3.450e+01 -27.730  < 2e-16 ***
## season8     -9.340e+02  3.451e+01 -27.061  < 2e-16 ***
## season9     -1.138e+03  3.452e+01 -32.972  < 2e-16 ***
## season10    -1.040e+03  3.453e+01 -30.122  < 2e-16 ***
## season11    -7.896e+02  3.497e+01 -22.577  < 2e-16 ***
## season12    -2.649e+02  3.498e+01  -7.571 9.72e-13 ***
## trend       -1.928e+00  4.479e-01  -4.304 2.51e-05 ***
## I(trend^2)   1.862e-02  1.676e-03  11.113  < 2e-16 ***
## s_break      6.060e+01  2.836e+01   2.137   0.0337 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109 on 223 degrees of freedom
## Multiple R-squared:  0.9423, Adjusted R-squared:  0.9387 
## F-statistic: 260.3 on 14 and 223 DF,  p-value: < 2.2e-16

Como puede observarse en el resumen del modelo anterior, la variable de ruptura estructural es estadísticamente significativa, con un nivel de 0,01. A la vez, en el caso de valores atípicos, el hot encoding puede aplicarse estableciendo una variable binaria igual a 1, siempre que un valor atípico se produzca y un valor de 0 de otro modo.

CONSEJO: Una vez que se haya practicado con un modelo de predicción con la función tslm y con el uso de variables externas, se tendrán que producir valores futuros de estas variables, ya que se van a utilizar como entrada del modelo de predicción.

Predicción de una serie con componentes de multiestacionalidad

Una de las principales ventajas del modelo de regresión, a diferencia de los modelos de series temporales tradicionales como ARIMA o Holt-Winters, es que proporciona una amplia gama de opciones de personalización y permite modelar y pronosticar datos de series temporales complejas como series con multiestacionalidad.

En los siguientes ejemplos, se hará uso de la serie UKgrid para demostrar el enfoque de predicción de una serie multiestacional con un modelo de regresión lineal.

La serie UKgrid

La serie UKgrid representa la demanda nacional de electricidad en el Reino Unido, y está disponible en el paquete UKgrid. Esta serie representa una serie temporal de alta frecuencia de media hora. Se usará la función extract_grid del paquete UKgrid para definir la serie, las características principales (por ejemplo, formato de datos, variables, frecuencia, etc.). Con esta función de transformación es posible agregar la frecuencia de la serie de media hora a una frecuencia más baja, (por ejemplo, por hora, diario o mensual).

Como el objetivo aquí es pronosticar la demanda diaria en los próximos 365 días, se establecerá la serie en frecuencia diaria usando la estructura data.frame:

library(UKgrid)

UKdaily <- extract_grid(type = "data.frame",
                        columns = "ND",
                        aggregate = "daily")

Se hará uso de la función principal para revisar las variables de la serie, dando como resultado:

head(UKdaily)

##    TIMESTAMP      ND
## 1 2005-04-01 1920069
## 2 2005-04-02 1674699
## 3 2005-04-03 1631352
## 4 2005-04-04 1916693
## 5 2005-04-05 1952082
## 6 2005-04-06 1964584

Como puede verse, esta serie tiene dos variables:

TIMESTAMP: un objeto de fecha utilizado como timestamp o índice de series

ND: La demanda neta de electricidad

También se utiliza la función ts_plot para trazar y revisar la estructura de la serie, como resultado la gráfica:

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

Como pudo verse anteriormente, la serie tiene una clara tendencia decreciente y un patrón de cadena estacional. Además, como se vio en el capítulo 6, Análisis de estacionalidad, esta serie tiene múltiples patrones de estacionalidad:

Diario: Un ciclo de 365 días al año

Día de la semana: Un ciclo de 7 días

Mensual: Efectuado por el tiempo

La evidencia de estos patrones se puede ver en el siguiente heatmap de la serie (desde 2016) usando la función ts_heatmap del paquete TSstudio, resultando la siguiente gráfica:

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

Como se puede ver en la serie de heatmap, la demanda global aumenta a lo largo de las semanas de la estación de invierno (por ejemplo, las semanas del calendario 1 a 12 y las semanas 44 a 52). Además, se puede observar el cambio de la serie durante los días de las semanas, ya que la demanda aumenta durante los días laborables 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 establece la serie como frecuencia una de diaria y se crean las dos siguientes características:

Indicador del día de la semana
Indicador del mes del año

Además, como es razonable asumir (esta suposición se confirmará con la función ACF una vez que se haya transformado la serie en un objeto ts) que la serie tiene una fuerte correlación con los retardos estacionales, se creará una variable de retardo de 365 observaciones. Además, se utilizará el paquete dplyr para crear las siguientes 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)

INFORMACIÓN: Es necesario recordar que el costo de usar una variable de retraso con una longitud de N es la pérdida de las primeras observaciones para N (ya que los retrasos de esas observaciones no se pueden generar a partir de la serie). Por lo tanto, se utilizan las funciones de filtro para eliminar las filas de la tabla que falta la variable lag365 (es decir, las primeras 365 observaciones).

Ahora, es necesario revisar la estructura de la tabla UKdaily después de añadir estas nuevas características, obteniendo la siguiente salida:

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 "Sun"<"Mon"<"Tue"<..: 7 1 2 3 4 5 6 7 1 2 ...
##  $ month    : Ord.factor w/ 12 levels "Jan"<"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 un formato ts (al menos para la serie), se convierte la serie en un objeto ts. Se utiliza 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:

start_date <- min(UKdaily$TIMESTAMP)
start <- c(year(start_date), yday(start_date))

Como pudo verse en el capítulo 3 The Time Series Object, se usa la función ts del paquete stats para establecer el objeto ts:

UK_ts <- ts(UKdaily$ND, 
            start = c(year(start_date), yday(start_date)),
            frequency = 365)

Después de transformar la serie en un objeto ts, se puede volver atrás y confirmar la suposición que se hizo sobre el nivel de correlación entre la serie y sus intervalos estacionales con la función ts_acf (una versión interactiva de la función acf del paquete TSstudio).

Se revisa la correlación de la serie con sus rezagos de los últimos cuatro años, obteniendo la gráfica de salida:

library(TSstudio)
data(USgas)
acf(UK_ts, lag.max = 365 * 4)

La trama ACF anterior confirma nuestra suposición, y la serie tiene una fuerte relación con los rezagos estacionales, en particular el rezago de 365, el primero.

CONSEJO: Como nota al margen, se puede estar seguro de que la serie también tiene una fuerte correlación con los retrasos semanales (es decir, retraso 7, 14, 21, y así sucesivamente). Sin embargo, en general, no se recomienda utilizar los retrasos que son más pequeños que el horizonte de pronóstico (por ejemplo, retraso 7, cuando el horizonte de previsión es 365), como se pueden prever esos rezagos, así como poder utilizarlos como insumo en el modelo. Esto implica un esfuerzo adicional, ya que tendrá que prever los retrasos. Además, se puede aumentar el sesgo de predicción utilizando valores previstos como entradas.

Ahora, después de haber creado las nuevas características para el modelo y establecer el objeto ts, está todo listo para dividir la serie de entrada y el objeto de características externas correspondientes que se crearon (UKdaily) en una partición de entrenamiento y prueba. Como el objetivo es pronosticar las siguientes 365 observaciones, y la longitud de la serie es lo suficientemente grande (2.540 observaciones), puede establecerse la partición de prueba a la longitud de las observaciones del horizonte de predicción-365.

Se establece h, una variable indicadora, a 365 y se utilizará para definir las particiones y, más tarde, el horizonte de predicción:

h <-  365

Como antes, se divide la serie en particiones de entrenamiento y pruebas con la función ts_split:

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

De manera similar, se dividen las características creadas para el modelo de regresión (las características estacionales y de retardo) en una partición de entrenamiento y prueba siguiendo el mismo orden exacto usado para el objeto ts correspondiente. Se utiliza la funcionalidad de índice data.frame para configurar la tabla UKdaily a particiones de entrenamiento y pruebas:

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

Entrenamiento y prueba del modelo de previsión

Una vez creadas las diferentes características del modelo, este se prepara para iniciar el proceso de entrenamiento del modelo de previsión. Se utilizará la partición de entrenamiento y se empezarán a entrenar los tres modelos siguientes:

Modelo de referencia: Es la regresión de la serie con el componente estacional y de tendencia en donde se utilizan las características incorporadas de la función tslm. Como se fija la frecuencia de la serie en 365, la característica estacional de la serie hace referencia a la estacionalidad diaria.
Modelo multiestacional: En este modelo se añaden los indicadores referentes a día de la semana y mes del año que permiten capturar la multiestacionalidad de la serie.
Un modelo multiestacional con un desfase estacional: Utilizando, además de los indicadores estacionales, la variable de retardo estacional.

La comparación de estos tres modelos se basará en los siguientes criterios:

Rendimiento del modelo en el conjunto de entrenamiento y prueba utilizando la puntuación MAPE
Visualización de los valores ajustados y previstos frente a los valores reales de la serie utilizando la función test_forecast

CONSEJO: Empezar con un modelo simplista (modelo de referencia) permitirá observar si el añadir nuevas características contribuye al rendimiento del modelo o si este debe ser evitado, esto debido a que añadir más características o complejidad al modelo no siempre produce mejores resultados.

Se empezará con el modelo de referencia, en donde se hará una regresión de la serie con sus componentes estacionales y de tendencia como se observa a continuación:

md_tslm1 <- tslm(train_ts ~ season + trend)

A continuación, se utilizará el modelo entrenado, md_tslm1, para predecir 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:

fc_tslm1 <- forecast(md_tslm1, h = h)

Para comparar el rendimiento del modelo en los conjuntos de entrenamiento y de prueba utilizando la función test_forecast se hace:

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

Se puede observar en el gráfico de rendimiento anterior que el modelo de referencia hace un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año. Por el contrario, no logra captar la oscilación relacionada con el día de la semana. La puntuación MAPE del modelo, como se observa en la salida de la siguiente función accuracy, es del 6,09% y del 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.071899

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

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

Dado que ahora se usan características de una fuente de datos externa, se deben especificar los datos de entrada con el argumento data. Ahora se repetirá el mismo proceso realizado anteriormente, al utilizar un pronóstico con el modelo entrenado y evaluar el rendimiento del modelo:

fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
              forecast.obj = fc_tslm2,
              test = test_ts)

Al observar el gráfico de rendimiento precedente del segundo modelo, se puede apreciar la contribución de las características estacionales en la previsión, puesto que el segundo modelo fue capaz de captar tanto la tendencia como los patrones multiestacionales de la serie. Esto también puede observarse en la puntuación del MAPE del modelo, que se redujo al 2,87% y al 5,23% en las particiones de entrenamiento y prueba, respectivamente:

accuracy(fc_tslm2, test_ts)

##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -9.812170e-13 70245.98 52146.79 -0.1738708 3.167605 0.4370853
## Test set     -1.764563e+04 80711.71 65373.21 -1.3715505 4.682071 0.5479470
##                   ACF1 Theil's U
## Training set 0.7513664        NA
## Test set     0.6075598   0.68445

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

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

El rendimiento del tercer modelo puede verse en el siguiente gráfico:

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

Es difícil observar en el gráfico de rendimiento del tercer modelo si hay una diferencia significativa con respecto al segundo. Por lo tanto, se revisa la puntuación MAPE del tercer modelo, donde se obtiene el siguiente resultado:

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

Los resultados del tercer modelo revelan una ligera mejora en la precisión del modelo, con un 2,81% en el conjunto de entrenamiento y un 5,07% en el conjunto de pruebas.

Model selection

En este punto debemos seleccionar un modelo. En este momento, está demostrado que el segundo y el tercer modelo obtienen mejores resultados que el primer modelo. Puesto que tanto el segundo como el tercer modelo obtuvieron una puntuación MAPE bastante similar, con una pequeña ventaja para el tercer modelo, conviene preguntarse si una mejora de menos del 0,2% en la tasa de error del conjunto de pruebas merece el coste de 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).

Depende.

No hay una respuesta directa a esta pregunta. Además, depende del objetivo de la previsión. Se recomienda tener en cuenta la siguiente prueba:

La primera pregunta a plantear en este caso es ¿Es la variable de retraso estadísticamente significativa? Si la variable no puede ser considerada estadísticamente significativa, no tiene sentido seguir con el debate y es mejor descartar la variable. En el caso del tercer modelo, podemos utilizar la función de resumen para observar el nivel de significación de la variable lag365 donde obtenemos:
```
summary(md_tslm3)$coefficients %>% tail(1)
```
```
##           Estimate Std. Error   t value     Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05
```
El valor p de la variable lag365 indica que la variable es estadísticamente significativa a un nivel de 0,001.

Del mismo modo, se puede aplicar una sola prueba ANOVA con la función anova del paquete stats, y comprobar si la variable adicional es significativa, obteniendo de esa manera:

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 ' ' 1

La prueba ANOVA indica de igual manera que la variable lag365 es estadísticamente significativa.

Prueba de respaldo: Se podría dar el caso de que el tercer modelo sea más acertado simplemente por casualidad y no debido a que la variable de retardo contribuya a la precisión del modelo, dado que la diferencia es relativamente pequeña. Por tanto, la prueba retrospectiva o de respaldo de ambos modelos podría ser útil para validar si la contribución de la variable de retardo es coherente a lo largo de varios periodos de prueba.

Información: Es posible aplicar métodos más sólidos para la selección de características, tales como la regresión escalonada, de cresta o de lazo. Aunque estos métodos no entran en el ámbito del libro, se recomienda leer sobre ellos. En el capítulo 12, Predicción con modelos de aprendizaje automático, se explorarán enfoques de regresión avanzados con modelos de aprendizaje automático, que incluyen métodos de selección de características.

En pro de la simplicidad, se usará el criterio de precisión y se seleccionará el tercer modelo para pronosticar las series. El próximo punto será volver a entrenar el modelo en todas las series y pronosticar los próximos 365 días de la siguiente manera:
```
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365, 
                 data = UKdaily)
```

Análisis de los residuos

Justo antes de finalizar la previsión, se realiza un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals, donde se obtiene:

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

Como se puede observar en el gráfico de resumen de residuos anterior, los residuales no son ruido blanco, pues existe cierta autocorrelación entre las series residuales y sus retardos. Esto es, técnicamente, una indicación de que el modelo no capturó todos los patrones o la información que existe en las series. Una manera de solucionar este problema es identificando variables adicionales que puedan representar la variación de los residuos. El principal obstáculo de este enfoque es que es complicado determinar las variables externas que pueden explicar la variación de los residuos y que también son viables a la hora de hacer previsiones. Por ejemplo, es lógico asumir que los modelos meteorológicos afectan a la demanda de electricidad, pero es difícil predecir esos modelos meteorológicos con un año de antelación.

Un enfoque alternativo el cual puede ser tenido en cuenta cuando los patrones sobran en los residuos del modelo, consiste en tratar los residuos del modelo como una serie de datos de tiempo separada y modelarla con otro modelo de previsión de series de tiempo. Un enfoque común es una regresión con error ARIMA, que se presentará en el capítulo 11, Previsión con modelos ARIMA.

Finalización de la previsión

Para finalizar el proceso, se utilizará el modelo entrenado seleccionado para pronosticar las futuras 365 observaciones. Como se utilizaron variables externas con la función tslm, será necesario generar sus valores futuros. Esto es algo relativamente sencillo, puesto que se utilizan variables deterministas. Por lo tanto, se creará un objeto data.frame con los valores de wday, month, y lag365 para las próximas 365 observaciones futuras. Un enfoque simplista consiste en crear primero las fechas correspondientes de las observaciones pronosticadas, y luego extraer de este objeto el día de la semana y el mes del año, como se observará a continuación:

UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1), 
                                       by = "day", 
                                       length.out = h))

A continuación, el paquete lubridate permite utilizar la variable date para crear las variables wday y month de la siguiente manera:

UK_fc_df$wday <- factor(lubridate::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)

Si se quiere crear una previsión, hay que utilizar la función forecast como se muestra a continuación:

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

Por último, pero no menos importante, se traza la previsión final con la función plot_forecast del paquete TSstudio, donde se obtiene:

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

Resumen

En este capítulo, se han presentado las aplicaciones de previsión del modelo de regresión lineal. Si bien el modelo de regresión lineal no fue ideado para manejar datos de series temporales, con una sencilla manipulación de características se puede transformar un problema de previsión en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal con respecto a otros modelos tradicionales de series temporales es la facilidad para incorporar variables y factores externos. variables y factores externos. No obstante, este modelo puede manejar series temporales con patrones de multiestacionalidad, como vimos con la previsión de la demanda de electricidad en el Reino Unido.

Parte 2

En esta parte se evaluarán los códigos anteriormente estudiados aplicados a una base de datos que proporciona información sobre la ocurrencia de sismos entre los años de 1965 a 2016.

terremotos<-read.csv("database.csv")
class(terremotos)

## [1] "data.frame"

library(dplyr)
terremotos <-terremotos  %>%
  arrange(Date)

terremotos <- ts(terremotos$Magnitude, start = c(1965,1), end = c(2016,12), frequency = 12)

class(terremotos)

## [1] "ts"

library(TSstudio)

data("terremotos")

## Warning in data("terremotos"): data set 'terremotos' not found

ts_plot(terremotos,
        title = "Terremotos registrados en Estados Unidos entre el año 1965 hasta el año 2016",
        Ytitle = "Magnitud",
        Xtitle = "Año")

En la serie de tiempo ejecutada anteriormente se observan los terremotos ocurridos entre el año de 1965 hasta 2016. Estos datos no presentan alguna frecuencia. Es decir, existe cierta aleatoriedad de sucesos.

ts_info(terremotos)

##  The terremotos series is a ts object with 1 variable and 624 observations
##  Frequency: 12 
##  Start time: 1965 1 
##  End time: 2016 12

Con el objetivo de descomponer la serie de tiempo y observar la tendencia y demás variables de la misma se utiliza:

ts_decompose(terremotos)

Se puede ver en el gráfico anterior que la serie no cuenta con una tendencia marcada a lo largo de los años. Por lo tanto, se puede afirmar que para este conjunto de datos hay aleatoriedad.

terremotos_df <- ts_to_prophet(terremotos)

La función transforma el objeto ts dentro de dos columnas de data.frame, en la cual, una de las columnas representan el tiempo y la otra columna los componentes de la serie:

head(terremotos_df)

##           ds   y
## 1 1965-01-01 6.5
## 2 1965-02-01 5.9
## 3 1965-03-01 5.6
## 4 1965-04-01 5.6
## 5 1965-05-01 6.0
## 6 1965-06-01 6.8

terremotos_df$trend <- 1:nrow(terremotos_df)

La regresión de la serie con su índice proporciona una estimación del crecimiento marginal de un mes a otro, ya que el índice se encuentra en un orden cronológico con un incremento constante.

La segunda característica de la regresión es el componente estacional. Debido a que se quiere medir la contribución de cada unidad de frecuencia a la oscilación de serie, se hace uso de una variable categórica para cada unidad de frecuencia. En el caso de la serie terremotos, las unidades de frecuencia representan los meses del año, por lo cual se crea una variable categórica con 12 categorías, en donde cada categoría es un mes del año. Por otro lado, se usa la función de month del paquete Lubridate para extraer los meses del año de la variable de fecha ds.

library(lubridate)

terremotos_df$seasonal <- month(terremotos_df$ds, label = T)

head(terremotos_df)

##           ds   y trend seasonal
## 1 1965-01-01 6.5     1      Jan
## 2 1965-02-01 5.9     2      Feb
## 3 1965-03-01 5.6     3      Mar
## 4 1965-04-01 5.6     4      Apr
## 5 1965-05-01 6.0     5      May
## 6 1965-06-01 6.8     6      Jun

h <- 12 

train <- terremotos_df[1:(nrow(terremotos_df) - h), ]

test <- terremotos_df[(nrow(terremotos_df) - h + 1):nrow(terremotos_df), ]

Después de crear los marcos de datos de entrenamiento y prueba, se puede revisar como el modelo de regresión toma cada uno de los componentes juntos y separados.

Modelando la tendencia de la serie y estacional componentes

En primer lugar, se modela la tendencia de la serie, en la cual, se hace una regresión de la serie con la variable de tendencia, en la partición de entrenamiento:

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

Se hace uso de la función summary para revisar los detalles del modelo:

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

Se prosigue a utilizar la función de forecast del paquete de estadísticos, el cual predice los valores de datos de entrada basados en un modelo determinado.

De igual manera, se usa md_trend para predecir los valores ajustados y previstos del modelo de tendencia anterior:

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

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

Se crea una función de utilidad que traza la serie y la salida del modelo, usando el paquete ploty:

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 = ""),
           yaxis = list(title = "Magnitud"),
           legend = list(x = 0.05, y = 0.95))
  return(p)
}

La función tiene los siguientes argumentos:

$\cdot$ data: los datos de entrada, un objeto data.frame que sigue la misma estructura que la de terremotos (incluyendo la variable yhat ).

$\cdot$ train: Conjunto de entrenamiento correspondiente que se usó para entrenar el modelo.

$\cdot$ test: Conjunto de pruebas correspondiente que se utilizó para evaluar el modelo pronóstico.

$\cdot$ tittle: título de la trama, se utiliza por defecto NULL.

Se establece las entradas de la función plot_lm con el modelo de salida:

plot_lm(data = terremotos_df, 
        train = train, 
        test = test,
        title = "Predicción del componente de tendencia de la serie")

En general, el modelo fué capaz de capturar el movimiento general de tendencia. A pesar de que la tendencia se mantiene uniforme a través de los años, puede notarse que la linea decrece muy lentamente. Para el análisis de comparación, se desea medir la tasa de error del modelo de entrenamiento y conjuntos de prueba: La norma se mide por medio de:

$MAPE=(\frac{1}{M})\sum|\frac{A_i-F_i}{A_i}|$

mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
                mean(abs(test$y - test$yhat) / test$y))

"Se obtiene la siguiente salida: "

## [1] "Se obtiene la siguiente salida: "

mape_trend

## [1] 0.05435221 0.04106426

El proceso de modelar y pronosticar el componente estacional sigue el mismo proceso de tendencia, al realizar una regresión de la serie con la variable estacional creada anteriormente:

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

Se prosigue a actualizar 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)

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

Como se observa en el gráfico anterior, el modelo falla a la hora de capturar adecuadamente la estructura del patrón estacional de la serie. Además, se puede denotar que hace falta la tendencia de la serie. Por lo tanto, antes de agregar la tendencia y los componentes estacionales, se estudia 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.05383805 0.04855883

Se puede decir que existe una alta tasa de error en el conjunto de pruebas relacionada con el componente de tendencia que no se incluye en el modelo anterior. Por otro lado, el siguiente paso a realizar es unir los dos componentes en un modelo y posterior a esto, pronosticar los valores característicos de la serie:

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 = terremotos_df, 
        train = train, 
        test = test,
        title = "Predicción del componente de tendencia de la serie")

Se mide 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.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 = terremotos_df, 
        train = train, 
        test = test,
        title = "Predicción del componente de tendencia de la serie")

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

La función tslm

Ahora se repetirá el ejemplo anterior y se pronosticarán las últimas 12 observaciones de la serie terremotos con la función tslm utilizando la tendencia, el cuadrado de la tendencia, y los componentes estacionales. En primer lugar, se divide la serie en particiones de entrenamiento y pruebas usando la función ts_split:

terremotos_split <- ts_split(terremotos, sample.out = h)

train.ts <- terremotos_split$train

test.ts <- terremotos_split$test

A continuación, se aplicará la misma fórmula que se usó para crear el anterior modelo de predicción md2 usando la función tslm:

library(forecast)

md3 <- tslm(train.ts ~ season + trend + I(trend^2))

Ahora se va a revisar el md3, el resultado de la función tslm, y compararla con el resultado de md2, obteniendo el resultado:

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

Eficiencia: no es necesario transformar la serie en un objeto data.frame y, además, ingeniería de características.
El objeto de salida soporta toda la funcionalidad del modelo de predicción (como las funciones de exactitud y checkresiduals) y paquetes de TSstudio (como las funciones test_forecast y plot_forecast)

Modelado de eventos individuales y eventos no estacionales

En algunas ocasiones, los datos de series temporales pueden o no contener patrones inusuales a través del tiempo. Los siguientes son ejemplos de tales eventos:

Valores atípicos: Evento(s) que están fuera de los patrones normales de la serie.

Ruptura estructural: es algún evento significativo que genera un cambio en los patrones históricos de la serie. Un ejemplo bastante común es un cambio en el crecimiento de la serie.
Eventos recurrentes no estacionales: sucede cuando un evento se repite de ciclo a ciclo, pero el momento en el que ocurren dichos cambios cambia de ciclo a ciclo. Un ejemplo común de este evento son las vacaciones de Pascua, las cuales ocurren de Marzo a Abril todos los años.

Por ejemplo, es posible observar en la trama de descomposición de la serie terremotos que fue expuesta anteriormente, que la tendencia de la serie tuvo una ruptura estructural alrededor del año 2010. Mientras que el crecimiento antes del 2010 fue relativamente nula, la pendiente de la tendencia cambió después, demostrando un crecimiento positivo. En este caso, es posible hacer uso de una variable binaria que sea igual a cero para las observaciones del año 2010 y un año después.

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

Ahora se usará la nueva característica para remodelar la serie terremotos:

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

Se usará la función de resumen para revisar la salida del modelo:

summary(md3)

## 
## Call:
## tslm(formula = terremotos ~ season + trend + I(trend^2) + s_break, 
##     data = terremotos_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

Predicción de una serie con componentes de multiestacionalidad

head(terremotos_df)

##           ds   y trend seasonal s_break
## 1 1965-01-01 6.5     1      Jan       0
## 2 1965-02-01 5.9     2      Feb       0
## 3 1965-03-01 5.6     3      Mar       0
## 4 1965-04-01 5.6     4      Apr       0
## 5 1965-05-01 6.0     5      May       0
## 6 1965-06-01 6.8     6      Jun       0

Como puede verse, esta serie tiene 4 variables:

ds: un objeto de fecha utilizado como timestamp o índice de series

Y: La magnitud neta de los terremotos

trend:

seasonal: El mes asociado a cada magnitud .

Usando la función ts_heatmap del paquete TSstudio, resultando la siguiente gráfica heatmap de la serie (desde 2016) :

library(TSstudio)
ts_heatmap (terremotos)

En donde podemos observar patrones de magnitud anuales distribuidos en los 12 meses del año .Es asi como para el mes de Agosto del año 1967 se registra el mayor evento de magnitud de la serie con magnitud de 7,9. Aunque no se puede establecer una tendencia para la magnitud de los terremmotosvemos como los eventos de mayor magnitud sísmica se concentran en marzo a Octubre de los años 1965-1970, de la misma manera como de 1993 a 2000 se presentan eventos de muy baja magnitud sísmica.

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

Para capturar los componentes estacionales de la serie, se establece la serie como frecuencia una de diaria y se crean las dos siguientes características:

Indicador del día de la semana
Indicador del mes del año

library(dplyr)
terremotos_dz<- terremotos_df%>%
  mutate(wday = wday(ds, label = TRUE),
         month = month(ds, label = TRUE),
         lag365 = dplyr::lag(trend,12
                             ))%>%
  filter(!is.na(lag365)) %>%
  arrange(ds)

Ahora, es necesario revisar la estructura de la tabla terremotos después de añadir estas nuevas características, obteniendo la siguiente salida:

str(terremotos_dz)

## 'data.frame':    612 obs. of  8 variables:
##  $ ds      : Date, format: "1966-01-01" "1966-02-01" ...
##  $ y       : num  5.5 5.7 6.1 6.3 6.7 5.5 6.7 6.3 5.5 7.2 ...
##  $ trend   : int  13 14 15 16 17 18 19 20 21 22 ...
##  $ seasonal: Ord.factor w/ 12 levels "Jan"<"Feb"<"Mar"<..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ s_break : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ wday    : Ord.factor w/ 7 levels "Sun"<"Mon"<"Tue"<..: 7 3 3 6 1 4 6 2 5 7 ...
##  $ month   : Ord.factor w/ 12 levels "Jan"<"Feb"<"Mar"<..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ lag365  : int  1 2 3 4 5 6 7 8 9 10 ...

start_date <- min(terremotos_dz$ds)

start <- c(year(start_date), yday(start_date))

Como pudo verse en el capítulo 3 The Time Series Object, se usa la función ts del paquete stats para establecer el objeto ts:

library(TSstudio)
terremotos_ts <- ts(terremotos_dz$y, 
            start = c(year(start_date), yday(start_date)),
            frequency = 12
            )

Se revisa la correlación de la serie con sus rezagos de los últimos cuatro años, obteniendo la gráfica de salida:

library(TSstudio)
data(terremotos_ts)

## Warning in data(terremotos_ts): data set 'terremotos_ts' not found

acf(terremotos_ts, lag.max = 12
    * 4)

La trama ACF anterior confirma nuestra suposición, y la serie tiene una fuerte relación con los rezagos estacionales, en particular el rezago de 365, el primero.

Ahora, después de haber creado las nuevas características para el modelo y establecer el objeto ts, está todo listo para dividir la serie de entrada y el objeto de características externas correspondientes que se crearon (terremotos) en una partición de entrenamiento y prueba. Como el objetivo es pronosticar las siguientes 365 observaciones, y la longitud de la serie es lo suficientemente grande (2.540 observaciones), puede establecerse la partición de prueba a la longitud de las observaciones del horizonte de predicción-365.

Se establece h, una variable indicadora, a 12 y se utilizará para definir las particiones y, más tarde, el horizonte de predicción:

h <-  12

Como antes, se divide la serie en particiones de entrenamiento y pruebas con la función ts_split:

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

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

Entrenamiento y prueba del modelo de previsión

Se utilizará la partición de entrenamiento y se empezarán a entrenar los tres modelos siguientes:

Modelo de referencia: Es la regresión de la serie con el componente estacional y de tendencia en donde se utilizan las características incorporadas de la función tslm. Como se fija la frecuencia de la serie en 365, la característica estacional de la serie hace referencia a la estacionalidad diaria.
Modelo multiestacional: En este modelo se añaden los indicadores referentes a día de la semana y mes del año que permiten capturar la multiestacionalidad de la serie.
Un modelo multiestacional con un desfase estacional: Utilizando, además de los indicadores estacionales, la variable de retardo estacional.

Se empezará con el modelo de referencia, en donde se hará una regresión de la serie con sus componentes estacionales y de tendencia como se observa a continuación:

library(TSstudio)
md_tslm1 <- tslm(train_ts ~ season + trend)

A continuación, se utilizará el modelo entrenado, md_tslm1, para predecir los próximos 12 días de la serie correspondientes a las observaciones de la partición de prueba, utilizando la función de previsión:

fc_tslm1 <- forecast(md_tslm1, h =h)

Para comparar el rendimiento del modelo en los conjuntos de entrenamiento y de prueba utilizando la función test_forecast se hace:

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

accuracy(fc_tslm1, test_ts)

##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set  0.0000000 0.4379175 0.3294690 -0.4959362 5.398051 0.7681513
## Test set     -0.1261184 0.2790590 0.2455092 -2.3655824 4.312106 0.5724005
##                      ACF1 Theil's U
## Training set -0.008509968        NA
## Test set     -0.553134413 0.7223825

Análisis de los residuos

Justo antes de finalizar la previsión, se realiza un análisis de los residuos del modelo seleccionado utilizando la función checkresiduals, donde se obtiene:

checkresiduals(terremotos)

## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Como se puede ver en el gráfico de resumen de residuos anterior, los residuos no son ruido blanco, ya que existe cierta autocorrelación entre las series de residuos y sus rezagos. Esto es técnicamente una indicación de que el modelo no capturó todos los patrones o la información que existe en las series. Una forma de abordar este problema es identificar variables adicionales que puedan explicar la variación de los residuos. El principal problema de este enfoque es que es difícil identificar variables externas que puedan explicar la variación de los residuos y que sean y que sean factibles de pronosticar. Por ejemplo, es difícil predecir patrones sismicos con un mes de antelación. Un enfoque alternativo, cuando los patrones sobran en los residuos del modelo, es tratar los de un modelo es tratar los residuos del modelo como datos de una serie temporal independiente y modelizarlos con otro modelo de previsión de series temporales.

Tercer informe: Análisis estádístico

Hernández Jhon, Llanos Guisell, Perdomo Juan, Rodriguez Laura

11/1/2021

Parte 1

Requisitos técnicos

Regresión lineal

Estimación de los coeficientes con el método OLS

Supuestos del OLS

Pronóstico con regresión lineal

Pronosticar la tendencia y los componentes estacionales

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

Modelando la tendencia de la serie y estacional componentes

La función tslm

Modelado de eventos individuales y eventos no estacionales

Predicción de una serie con componentes de multiestacionalidad

La serie UKgrid

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

Entrenamiento y prueba del modelo de previsión

Model selection

Análisis de los residuos

Finalización de la previsión

Resumen

Parte 2

Modelando la tendencia de la serie y estacional componentes

La función tslm

Modelado de eventos individuales y eventos no estacionales

Predicción de una serie con componentes de multiestacionalidad

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

Entrenamiento y prueba del modelo de previsión

Análisis de los residuos