PARTE 1

Características de la ingeniería de los componentes de la serie

Anterior a generar las entradas de regresión que representan la tendencia de la serie y los elementos estacionales. Primero tenemos que comprender su composición. En los próximos ejemplos el proceso de construcción de novedosas propiedades desde series existentes usando la serie USgas. Carguemos otra vez la serie a partir del paquete TS studio y representémosla con la funcion ts_plot:

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

Además, vamos a revisar las principales características de la serie utilizando 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 ver en el gráfico de la serie, y como vimos en los capítulos anteriores, USgas es una serie mensual con un intenso elemento estacional mensual y una línea de tendencia bastante estable. Tenemos la posibilidad de explorar la composición de los elementos de la serie con la funcionalidad ts_decompose:

ts_decompose(USgas)

Se puede ver en el gráfico anterior que la tendencia de la serie es bastante plana entre 2000 y 2010, y tiene un aumento bastante lineal hacia adelante. Por consiguiente, la tendencia general entre 2000 y 2018 no es estrictamente lineal. Este es un dato fundamental que nos permitirá conceptualizar el ingreso de la tendencia para el modelo de regresión.

Anterior a usar la funcionalidad lm, la capacidad de regresión lineal incorporada en R del paquete stats tendremos que cambiar la serie de un objeto ts a un objeto data.frame.Por consiguiente, usaremos la funcionalidad ts_to_prophet del paquete TSstudio:

USgas_df <- ts_to_prophet(USgas)

La funcionalidad transforma el objeto ts en 2 columnas de data.frame, donde ambas columnas representan los elementos temporales y numéricos de la serie, respectivamente:

head(USgas_df)

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

Luego de cambiar las series en un objeto data.frame, tenemos la posibilidad de comenzar a generar las propiedades de acceso de la regresión. La primera característica que crearemos es la tendencia de la serie. Un enfoque basico para edificar la variable de tendencia es indexar las visualizaciones de la serie en orden cronológico:

USgas_df$trend <- 1:nrow(USgas_df)

La regresión de la serie con el índice de la serie da una estimación del incremento marginal de mes a mes, debido a que el índice está en orden cronológico con incrementos constantes. La segunda característica que deseamos producir es el elemento estacional. Ya que deseamos medir la contribución de cada unidad de frecuencia a la oscilación de la serie, usaremos una variable categórica para cada unidad de frecuencia. En la situación de la serie USgas, las unidades de frecuencia representan los meses del año y, por consiguiente, crearemos una variable categórica con 12 categorías, todas las cuales corresponde a un mes específico del año. Usaremos la funcionalidad month del paquete lubridate para sustraer el mes del año de la variable de fecha ds:

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

USgas_df$seasonal <- factor(month(USgas_df$ds, label = T), ordered = FALSE)

Usamos la funcionalidad factor para transformar la salida de la funcionalidad month en una variable categórica no ordenada. Comprobemos ahora el marco de datos luego de aumentar las novedosas funcionalidades:

head(USgas_df)

##           ds      y trend seasonal
## 1 2000-01-01 2510.5     1      ene
## 2 2000-02-01 2330.7     2      feb
## 3 2000-03-01 2050.6     3      mar
## 4 2000-04-01 1783.3     4      abr
## 5 2000-05-01 1632.9     5      may
## 6 2000-06-01 1513.1     6      jun

Finalmente, sin embargo, no menos fundamental, previo a comenzar a hacer la regresión de la serie con aquellas propiedades, la dividiremos en una partición de entrenamiento y otra de prueba. Estableceremos los últimos 12 meses de la serie como partición de prueba:

h <- 12 # setting a testing partition length
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]

Ahora, luego de producir los marcos de datos de entrenamiento y de prueba, vamos a comprobar cómo el modelo de regresión capta todos los elementos por separado y en grupo.

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

Primero modelaremos la tendencia de la serie realizando una regresión de la serie con la variable de tendencia, sobre la partición de entrenamiento:

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

Utilizaremos la función summary para comprobar 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 puede verse en el resultado de la regresión anterior, el coeficiente de la variable de tendencia es estadísticamente importante a un grado de 0,001. No obstante, el R-cuadrado ajustado de la regresión es bastante baja, lo cual generalmente tiene sentido, debido a que la mayoría de la alteración de la serie está relacionada con el jefe estacional, como hemos observado en los gráficos anteriores.

Como se puede mirar en el resultado de la regresión anterior, la cuarta columna representa el costo p de todos los coeficientes del modelo. El costo p otorga la posibilidad de que rechacemos la premisa nula ya que es verdadera, o el error de tipo I. Por consiguiente, para el costo p menor que α, el costo umbral, rechazaremos la premisa nula con un grado de importancia α, donde los valores normales de α son 0,1, 0,05, 0,01, etc.

Como constantemente, se ofrece contextualizar los números con la visualización de los datos. visualización de los datos. Por consiguiente, usaremos el modelo que hemos diseñado para pronosticar los valores ajustados en la partición de entrenamiento y los valores pronosticados en la partición de prueba. La funcionalidad predict del paquete stats, como su nombre sugiere, predice los valores de un dato de ingreso basándose en un modelo dado.

La usaremos para profetizar los valores ajustados y previstos del modelo de tendencia que hemos entrenado antes:

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

Vamos a crear una función de utilidad que traza las series y la salida del modelo, utilizando el paquete plotly:

library(plotly)

## Warning: package 'plotly' was built under R version 4.1.1

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     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:

-data: Los datos de entrada, un objeto data.frame que sigue la misma estructura que la de USgas_df (incluyendo la variable yhat). -train: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo. -test: Igualmente, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de previsión. -title: El título de la parcela, por defecto, establecido en NULL.

Establezcamos las entradas de la función plot_lm con la salida del modelo:

plot_lm(data = USgas_df,

train = train,

test = test,

title = "Predicting the Trend Component of the
Series")

El modelo pudo capturar el movimiento general, aunque es posible que una tendencia lineal no capture la ruptura estructural en la tendencia alrededor de 2010.

Finalmente, para el análisis comparativo, queremos medir la tasa de error del modelo tanto en el entrenamiento como en el conjunto de prueba:

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

El proceso de modelado y pronóstico del componente estacional sigue el mismo proceso que aplicamos con la tendencia, realizando una regresión de la serie con la variable estacional que creamos anteriormente:

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

Repasemos los detalles del modelo:

summary(md_seasonal)

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

Debido a que devolvemos la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías que están incluidas en los valores de pendiente. Como puede verse en el resumen de la regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. También se puede observar que el R-cuadrado ajustado del modelo estacional es ligeramente superior al del modelo de tendencia (0,78 versus 0,1).

Antes de trazar el modelo ajustado y los valores de previsión con la función plot_lm, vamos a actualizar los valores de yhat con la función predict:

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

Ahora podemos utilizar la función plot_lm para visualizar el modelo ajustado y los valores previstos:

plot_lm(data = USgas_df,

train = train,

test = test,

title = "Predicting the Seasonal Component of the
Series")

Como puede verse en el gráfico anterior, el modelo hace un buen trabajo al capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que falta la tendencia de la serie. Antes de agregar los componentes de tendencia y estacional, evaluemos el desempeño 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

La alta tasa de error en el conjunto de prueba está relacionada con el componente de tendencia que no se incluyó en el modelo. En el siguiente paso, los dos componentes se combinan en un modelo y se pronostican los valores característicos de la serie:

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

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

summary(md1)

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

El retorno de la serie con la tendencia y los componentes estacionales juntos proporciona un aumento adicional del R-cuadrado ajustado del modelo, que pasa de 0,78 a 0,91. Lo cual puede observarse en el gráfico de los resultados del modelo:

train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend and Seasonal Components of the
Series")

Midamos la puntuación MAPE del modelo en las particiones de entrenamiento y de 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 la tendencia y los componentes estacionales conduce a una mejora significativa tanto en la calidad del ajuste del modelo como en la precisión. Sin embargo, si se mira el gráfico del ajuste del modelo y el pronóstico, se puede ver que la tendencia del modelo es demasiado lineal y no muestra la ruptura estructural de la tendencia de la serie. Este es el punto en el que agregar un componente polinomial al modelo podría mejorar aún más la precisión del modelo.

Una técnica sencilla para capturar una tendencia no lineal es agregar un componente polinomial a la tendencia de la serie para capturar la curvatura de la tendencia a lo largo del tiempo. Usamos el argumento I que nos permite realizar operaciones matemáticas en cualquiera de los objetos de entrada. Entonces usamos este argumento para agregar un segundo grado del polinomio para la entrada de tendencia:

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

El resumen del modelo puede verse como sigue:

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

Agregar el polinomio de segundo grado al modelo de regresión no implicó una mejora significativa en la bondad de ajuste del modelo. En el otro modelo, como se ve en el siguiente gráfico de salida del modelo, este simple cambio en la estructura del modelo nos permite captar el desglose estructural de la tendencia a lo largo del tiempo:

train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend (Polynomial) and Seasonal Components
of the Series")

Como observamos en el modelo que sigue a la puntuación MAPE, la exactitud del modelo mejora significativamente al añadir la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas se redujo del 9,2% al 4,5%:

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

## [1] 0.03688770 0.04212618

La función tslm

Por el momento, hemos visto el proceso manual de transformación de un objeto ts a un formato de modelo de previsión de regresión lineal de regresión lineal. La función tslm del paquete forecast brinda una función integrada para transformar un objeto ts en un modelo de previsión de regresión lineal. Con la función tslm, puede establecer el componente de regresión junto con otras características.

Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie USgas con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. En primer lugar, vamos a dividir la serie en particiones de entrenamiento y de prueba utilizando la función ts_split.

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

A continuación, aplicaremos la misma fórmula que utilizamos para crear el modelo de previsión md2 anterior utilizando la función tslm:

library(forecast)

## 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 comparémosla con la salida de md2:

summary(md3)

## 
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.635e+03  3.224e+01  81.738  < 2e-16 ***
## season2     -3.004e+02  3.540e+01  -8.487 3.69e-15 ***
## season3     -4.841e+02  3.540e+01 -13.676  < 2e-16 ***
## season4     -9.128e+02  3.540e+01 -25.787  < 2e-16 ***
## season5     -1.099e+03  3.540e+01 -31.033  < 2e-16 ***
## season6     -1.118e+03  3.540e+01 -31.588  < 2e-16 ***
## season7     -9.547e+02  3.540e+01 -26.968  < 2e-16 ***
## season8     -9.322e+02  3.541e+01 -26.329  < 2e-16 ***
## season9     -1.136e+03  3.541e+01 -32.070  < 2e-16 ***
## season10    -1.046e+03  3.541e+01 -29.532  < 2e-16 ***
## season11    -8.001e+02  3.590e+01 -22.286  < 2e-16 ***
## season12    -2.618e+02  3.590e+01  -7.293 5.95e-12 ***
## trend       -1.270e+00  4.472e-01  -2.840  0.00494 ** 
## I(trend^2)   1.713e-02  1.908e-03   8.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9341 
## F-statistic: 246.1 on 13 and 212 DF,  p-value: < 2.2e-16

Como puede observar en la salida anterior, ambos modelos (md2 y md3) son idénticos.

El uso de la función tslm tiene varias ventajas, en lugar de establecer manualmente un modelo de regresión para las series con la función lm:

La eficiencia no requiere transformar las series a un objeto data.frame y ingeniería de características.
El objeto de salida soporta toda la funcionalidad de forecast (como las funciones accuracy y checkresiduals) y de los paquetes de TSstudio (como las funciones test_forecast y plot_forecast).

Modelización de eventos únicos y eventos no estacionales

En algunos casos, los datos de series de tiempo pueden contener o no patrones inusuales que se repiten con el tiempo. Los siguientes son ejemplos de estos tipos de eventos:

Valores atípicos: Un solo evento o eventos que están fuera del patrón normal de la serie.
Rotura estructural: hecho significativo que cambia los modelos históricos de la serie. Un ejemplo común es un cambio en el crecimiento de la serie.
Eventos no estacionales recurrentes: un evento que se repite de un ciclo a otro, pero su momento cambia de un ciclo a otro. Un ejemplo común de este tipo de eventos son las vacaciones de Semana Santa, que ocurren todos los años alrededor de marzo / abril.

Dado que este tipo de evento no se expresa en el modelo de regresión, sesgará los coeficientes estimados, ya que el modelo evaluará este tipo de evento junto con los eventos regulares de la serie. El uso de variables de código activo, binarias o de marca podría ayudar al modelo a ignorar este tipo de eventos o ajustar los coeficientes del modelo en consecuencia.

En el diagrama de descomposición de la serie de USgas de EE. UU. Que se muestra arriba, se puede ver, por ejemplo, que la tendencia de la serie alrededor de 2010 mostró una ruptura estructural. Si bien el crecimiento fue relativamente plano antes de 2010, la pendiente de la tendencia cambió posteriormente con un crecimiento positivo. En este caso, podemos usar una variable binaria que sea cero para las observaciones antes de 2010 y un año después.

La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejemplo muestra el proceso de creación de una variable binaria externa que es igual a 0 antes del año 2010 y a 1 después, utilizando 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 summary 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 se puede observar en el resumen del modelo anterior, la variable de rotura estructural es estadísticamente significativa con un nivel de 0.01. La codificación en caliente también se puede aplicar a valores atípicos o días festivos estableciendo una variable binaria igual a 1 siempre que sea un valor atípico o un evento no estacional recurrente 0.

De lo contrario, tenga en cuenta que después de haber entrenado un modelo de pronóstico con el La función tslm utiliza variables externas que deben generar valores futuros de estas variables, ya que se utilizarán como entrada para el pronóstico.

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

Una de las principales ventajas del modelo de regresión sobre los modelos tradicionales de series de tiempo como ARIMA o HoltWinters es que ofrece una amplia gama de opciones de personalización y nos permite compilar datos de series de tiempo complejas, como series con multiestacionalidad. .

En los siguientes ejemplos usamos la serie UKgrid para demostrar la predicción de una serie de varias estaciones con un modelo de regresión lineal.

La serie UKgrid

La serie UKgrid representa las necesidades de electricidad de la red nacional del Reino Unido y está disponible en el paquete UKgrid. Esta serie representa una serie temporal de datos de alta frecuencia con una frecuencia de media hora. Usamos la función extract_grid del paquete UKgrid para definir la serie, las características principales (por ejemplo, el formato de los datos, las variables de frecuencia, etc.). Esta función de transformación nos permite agregar la serie de media hora a una frecuencia más baja, por ejemplo horaria, diaria o mensual. Dado que nuestro objetivo es pronosticar la demanda diaria durante los próximos 365 días, utilizaremos la estructura data.frame para llevar la serie a la frecuencia diaria:

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:

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 ver, esta serie tiene dos variables:

TIMESTAMP: Un objeto de fecha utilizado como marca de tiempo o índice de la serie
ND: La demanda neta de electricidad

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

ts_plot(UKdaily,

title = "The UK National Demand for
Electricity",

Ytitle =
"MW",

Xtitle =
"Year")

Como puede observarse en el gráfico anterior, la serie tiene una clara tendencia a la baja y una cadena de patrones estacionales. Como vimos en el Capítulo 6, Análisis de estacionalidad, esta serie exhibe varios 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: Afectado por el clima

La evidencia de esos patrones puede verse en el siguiente mapa de calor de la serie desde 2016 utilizando la función ts_heatmap del paquete TSstudio:

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

Como se puede ver en el mapa de calor de la serie, la demanda global aumenta durante las semanas de invierno (por ejemplo, las semanas del calendario 1-12 y las semanas 44-52). Además, el cambio de las filas se puede observar los días laborables, ya que la demanda aumenta en los días laborables de la semana y disminuye los fines de semana.

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

Para captar los componentes estacionales de la serie, estableceremos la serie como diaria y crearemos las siguientes dos características:

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

También como es razonable suponer (hacemos esta suposición con la función ACF, una vez que transformamos la serie en un objeto ts) que la serie tiene una fuerte correlación con los rezagos estacionales, creamos una variable con un rezago de 365 observaciones. Usamos el paquete dplyr para crear estas 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)

Recuerde que el costo de usar una longitud de N variable de retardo es la pérdida de las primeras N observaciones (ya que los retrasos de estas observaciones no pueden estar fuera de secuencia). Por lo tanto, usamos las funciones de filtro para eliminar las filas de la tabla en las que falta la variable lag365 (es decir, las primeras 365 observaciones).

Ahora que hemos agregado estas nuevas características, echemos un vistazo a la estructura de la tabla UKdaily:

str(UKdaily)

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

Dado que la entrada de la función tslm tiene que estar en formato ts (al menos para la cadena), convertimos la cadena en un objeto ts. Usamos 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 inicio 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)

Después de convertir la serie en un objeto ts, podemos usar la función ts_acf (una versión interactiva de la función acf del paquete TSstudio) para confirmar nuevamente la suposición que hicimos sobre el grado de correlación entre la serie y sus rezagos estacionales. Comprobaremos la correlación de la serie con sus rezagos en los últimos cuatro años:

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

El gráfico ACF anterior confirma nuestra suposición, y la serie tiene una fuerte relación con los rezagos estacionales, en particular el rezago 365, el primer rezago.

Como nota al margen, puede estar seguro de que la serie también tiene una fuerte correlación con los retrasos semanales (es decir, retrasos 7, 14, 21, etc.). Sin embargo, en general, no se recomienda utilizar retrasos que sean menores que el horizonte de pronóstico (por ejemplo, retraso 7 si el horizonte de pronóstico es 365) porque necesita predecir estos retrasos para usarlos como entrada para el modelo de pronóstico. Esto requiere un esfuerzo adicional ya que necesita predecir estos retrasos. Además, puede aumentar la variación del pronóstico porque estamos utilizando valores de pronóstico como entradas.

Ahora que hemos creado las nuevas características para el modelo y establecido el objeto ts, podemos dividir la serie de entrada que creamos y el objeto de característica externa correspondiente (UKdaily) en una partición de entrenamiento y una partición de prueba. Dado que nuestro objetivo es predecir los próximos 365 y la longitud de la serie es lo suficientemente grande (2.540 observaciones), podemos permitirnos establecer la partición de prueba a la longitud del horizonte de pronóstico - 365 observaciones. Estableceremos h, una variable indicadora, en 365 y así definiremos las particiones y luego el horizonte de pronóstico:

h <- 365

Como antes, dividiremos la serie en una partición de entrenamiento y otra de prueba con la función ts_split que se utiliza para dividir la serie en una partición de entrenamiento y otra de prueba:

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

De manera similar, necesitamos dividir las características que construimos para el modelo de regresión (las características de temporada y de retraso) en una partición de entrenamiento y una partición de prueba, siguiendo exactamente el orden que usamos para el objeto ts correspondiente. Usamos la funcionalidad de índice data.frame para establecer la tabla UKdaily en las particiones de entrenamiento y prueba:

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

Entrenamiento y prueba del modelo de previsión

Una vez que se han creado las diferentes características del modelo, estamos listos para comenzar el proceso del modelo de pronóstico. Usaremos la partición de entrenamiento y comenzaremos a entrenar los siguientes tres modelos:

Modelo de referencia: regresión de componente estacional y tendencia usando la funcionalidad incorporada de la función tslm. Dado que la frecuencia de la serie es 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
Modelo multiestacional: agregue el día de la semana y el mes del año para capturar la multiestacionalidad de la serie.
Modelo multiestacional con retraso estacional: uso, además de indicadores estacionales, de la variable de retraso 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.
Mostrar valores ajustados y predichos versus valores reales de la serie usando la función test_forecast

Comenzar con un modelo simplista (modelo de referencia) nos permitirá ver si agregar nuevas características contribuye al desempeño del modelo o si debemos evitarlo, ya que agregar más funcionalidad o complejidad al modelo no siempre da mejores resultados.

Partamos del modelo de referencia, retrocediendo la serie con sus componentes estacionales y tendenciales:

md_tslm1 <- tslm(train_ts ~ season + trend)

A continuación, usamos el modelo entrenado md_tslm1 para predecir los próximos 365 días seguidos de acuerdo con las observaciones de la partición de prueba con forecast:

fc_tslm1 <- forecast(md_tslm1, h = h)

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

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

En el gráfico de rendimiento anterior, podemos ver que el índice de referencia hace un excelente trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año al capturar tanto la tendencia de la serie como la estacionalidad del día del año. . Por otro lado, no se registra la vibración relacionada con el día de la semana. Como podemos ver en el resultado de la siguiente función de accuracy, la puntuación MAPE del modelo es 6.09% y 7.77% en la partición de entrenamiento o prueba:

accuracy(fc_tslm1, test_ts)

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

Ahora intentaremos mejorar la precisión del modelo, añadiendo al modelo 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 + month,
data = train_df)

Como ahora estamos utilizando características de una fuente de data externa, tenemos que especificar los datos de entrada con el argumento data. Ahora repetiremos el mismo proceso que antes, utilizando una con el modelo entrenado y evaluando el rendimiento del modelo:

Si miramos el gráfico de desempeño anterior del segundo modelo, podemos observar la contribución de las características estacionales al pronóstico, ya que el segundo modelo fue capaz de capturar tanto la tendencia como los patrones multiestacionales de la serie. Esto también se puede ver en la puntuación MAPE del modelo, que se redujo a 2,87% y 5,23% en las particiones de entrenamiento y prueba:

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, pero no menos importante, añadamos la variable de retardo al modelo, y repitamos 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)

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 ver en el gráfico de rendimiento del tercer modelo si hay una diferencia significativa con el segundo modelo. A continuación, veamos la puntuación MAPE del tercer modelo:

accuracy(fc_tslm3, test_ts)

##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -1.175109e-11 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set     -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
##                   ACF1 Theil's U
## Training set 0.7500146        NA
## Test set     0.6094666 0.6925767

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

Selección del modelo

El momento de elegir un modelo es ahora. En este punto, está claro que el segundo y tercer modelo funcionan mejor que el primer modelo. Dado que tanto el segundo como el tercer modelo obtienen puntajes similares, con una pequeña ventaja para el tercer modelo, tenemos que preguntarnos si una mejora en la tasa de error del conjunto de prueba de menos del 0.2% reduciría el costo del valor al usar la prueba. variable set.lag (es decir, la pérdida de 365 observaciones y el costo agregado de un grado de libertad para el modelo).

Depende.

No hay una respuesta sencilla a esta pregunta. También depende del objetivo del pronóstico. Se recomienda considerar la siguiente prueba:

La primera pregunta en este caso es: ¿Es la variable de rezago estadísticamente significativa? Si la variable no es estadísticamente significativa, no tiene sentido continuar la discusión y es mejor abandonar la variable. En el caso del tercer modelo, podemos usar la función summary para observar el nivel de significancia de la variable lag365:

summary(md_tslm3)$coefficients %>% tail(1)

##           Estimate Std. Error   t value     Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-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, podemos aplicar una única prueba ANOVA con la función anova del paquete stats, y comprobar si la variable adicional es significativa:

anova(md_tslm3)

## Analysis of Variance Table
## 
## Response: train_ts
##             Df     Sum Sq    Mean Sq    F value    Pr(>F)    
## season     364 1.4279e+14 3.9227e+11    73.5348 < 2.2e-16 ***
## trend        1 7.2634e+13 7.2634e+13 13615.7078 < 2.2e-16 ***
## wday         6 4.5009e+13 7.5016e+12  1406.2214 < 2.2e-16 ***
## month       11 1.3721e+11 1.2473e+10     2.3382  0.007266 ** 
## lag365       1 8.1432e+10 8.1432e+10    15.2650  9.49e-05 ***
## Residuals 4190 2.2352e+13 5.3345e+09                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Backtesting: Es posible que el tercer modelo solo sea más preciso por casualidad y no porque la variable de rezago contribuya a la precisión del modelo, ya que la diferencia es relativamente pequeña. Por lo tanto, el backtesting de ambos modelos podría ayudar a validar si la contribución de las variables de rezago es consistente durante diferentes períodos de prueba. Te dejo el backtesting de ambos modelos como un ejercicio divertido.

Se pueden utilizar métodos más robustos para la selección de características, como B. Regresión por pasos, regresión de cresta o regresión de lazo. Aunque estos métodos están más allá del alcance de este libro, se recomienda que lo lea. En el Capítulo 12, Predicción mediante modelos de aprendizaje automático, examinaremos los enfoques de regresión mediante modelos de aprendizaje automático que incorporan métodos de selección de características.

Para mas simplicidad, nos atenemos al criterio de precisión y elegimos el tercer modelo para predecir la serie. El siguiente paso es volver a entrenar el modelo para todas las filas y pronosticar los próximos 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:

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 ver en el gráfico de resumen de residuos anterior, los residuos no son ruido blanco porque existe cierta autocorrelación entre la serie residual y sus retrasos. Esto es técnicamente una indicación de que el modelo no captura todos los patrones o la información que está en la serie. Una forma de abordar este problema es identificar variables adicionales que puedan explicar la variación en 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 puedan predecirse. Por ejemplo, es razonable suponer que los patrones climáticos afectan la demanda de electricidad, pero es difícil predecir estos patrones climáticos con un año de anticipación.

Si se dejan patrones en los residuos del modelo, un enfoque alternativo es tratar los residuos del modelo como datos de series de tiempo separados y modelarlos con un modelo de pronóstico de series de tiempo diferente. Modelo de pronóstico. Un enfoque común es la regresión de errores ARIMA, que presentamos en el Capítulo 11, Predicción con modelos ARIMA.

Finalización de la previsión

Terminemos el proceso y usemos el modelo entrenado que seleccionamos para predecir 365 observaciones futuras. Dado que usamos variables externas con la función tslm, necesitamos generar sus valores futuros. Esto es relativamente fácil porque usamos variables deterministas. Entonces creamos un data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras. Un enfoque simplificado es crear primero las observaciones pronosticadas y luego extraer el día de la semana y el mes del año de este objeto:

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

A continuación, podemos utilizar la variable 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)

Utilicemos la función forecast para crear la previsión:

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

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

plot_forecast(UKgrid_fc,

title = "The UK National Demand for Electricity
Forecast",

Ytitle =
"MW",

Xtitle =
"Year")

Resumen

En este capítulo, hemos presentado las aplicaciones de previsión del modelo de regresión lineal. Aunque el modelo de regresión lineal no se diseñó para manejar datos de series temporales, con con una simple ingeniería de características podemos transformar un problema de previsión en un problema de regresión lineal. en un problema de regresión lineal. La principal ventaja del modelo de regresión lineal con respecto a otros de series temporales tradicionales es la capacidad del modelo de incorporar variables y factores externos. variables y factores externos. No obstante, este modelo puede manejar series temporales con patrones de multiestacionalidad como hemos visto con la previsión de la demanda de electricidad en el Reino Unido.

PARTE 2

Lectura de datos geologicos

library(readr)
Datos_sismo <- read.csv("C:\\Users\\Jhona\\OneDrive\\Escritorio\\database.csv")
class(Datos_sismo)

## [1] "data.frame"

library(dplyr)
terremoto <-Datos_sismo  %>%
  arrange(Date)

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

class(Datos_sismo)

## [1] "ts"

Propiedades de la ingeniería de los componentes de la serie

Antes de generar elementos de regresión que representen una tendencia en elementos de matriz y estacionales y entender su composición a partir de lotes existentes utilizando usando la serie Datos_sismo.

Recarguemos la serie del paquete TS studio y representémosla con la funcionalidad ts plot.

library(TSstudio)

data(Datos_sismo)

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

ts_plot(Datos_sismo,

title =
"Magnitud de los Sismos desde 1965 a 2016",

Ytitle =
"Magnitud",

Xtitle =
"Year")

Igualmente,revisaremos características primordiales de la serie utilizando la función ts_info

ts_info(Datos_sismo)

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

En los capitulos anteriores Datos_sismo es una serie mensual con un intenso elemento estacional mensual y una línea de tendencia bastante estable, esto se puede evidenciar en como se ve en el gráfico de la serie,

Se tiene la posibilidad de inspeccionar la estructura de los elementos de la serie con la funcionalidad ts_decompose

ts_decompose(Datos_sismo)

La tendencia de la serie estacional aumenta a medida que aumenta el mes del año especificado hasta alcanzar su valor máximo, y comenzará a descender hasta principios del año siguiente, que se aplica a los 53 años presentados. La tendencia general entre 1965 y 2016 no es lineal.

Por eso,vamos a usar la funcionalidad ts_to_prophet del paquete TSstudio

Datos_sismo_df <- ts_to_prophet(Datos_sismo)

La funcionalidad transforma el objeto ts en 2 columnas de data.frame, donde ambas columnas representan los elementos temporales y numéricos de la serie, respectivamente. El código es el siguiente:

head(Datos_sismo_df)

##           ds   y
## 1 1965-01-01 6.0
## 2 1965-02-01 5.8
## 3 1965-03-01 6.2
## 4 1965-04-01 5.8
## 5 1965-05-01 5.8
## 6 1965-06-01 6.7

Después de modificar la serie en el objeto data.frame, es posible que comencemos a generar atributos de acceso de regresión. La primera característica que crearemos es la tendencia de la serie. Un método irremplazable o cambiante para construir variables de tendencia y visualización de series de índices en orden cronológico

Datos_sismo_df$trend <- 1:nrow(Datos_sismo_df)

La regresión de la secuencia al índice de secuencia da una estimación del segundo aumento cada mes, porque el índice se realiza en orden cronológico y continúa aumentando. La segunda característica que desea describir es la estacionalidad. Dado que queremos medir la contribución de cada unidad de frecuencia a la oscilación en serie, usaremos una variable categórica para cada unidad de frecuencia. En el caso de la serie Datos_sismo, la unidad de frecuencia representa el mes del año, por lo que crearemos una variable categórica con 12 categorías, correspondientes a un mes del año. Además, usaremos el paquete lubridate para restar el mes en el año de la variable de fecha ds

library(lubridate)
Datos_sismo_df$seasonal <- factor(month(Datos_sismo_df$ds, label = T), ordered = FALSE)

Ahora usaremos la funcionalidad factorial para variar la salida de la funcionalidad del mes en una variable categórica no ordenada. Se verifica el marco de datos luego de aumentar las funcionalidades

head(Datos_sismo_df)

##           ds   y trend seasonal
## 1 1965-01-01 6.0     1      ene
## 2 1965-02-01 5.8     2      feb
## 3 1965-03-01 6.2     3      mar
## 4 1965-04-01 5.8     4      abr
## 5 1965-05-01 5.8     5      may
## 6 1965-06-01 6.7     6      jun

Antes comenzar a hacer la regresión de la serie con dichas propiedades estableceremos los últimos 12 meses de la serie como partición de prueba:

h <- 12 # setting a testing partition length
train <- Datos_sismo_df[1:(nrow(Datos_sismo_df) - h), ]
test <- Datos_sismo_df[(nrow(Datos_sismo_df) - h + 1):nrow(Datos_sismo_df), ]

Luego de producir los marcos de datos de entrenamiento y de prueba,comprobramos si el modelo de regresión capta todos los elementos por separado y en grupo.

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

Primero configuraremos la tendencia de la serie haciendo una regresión de la serie con la variable de tendencia, sobre la partición de entrenamiento, el código resultante es:

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

Para el código utilizaremos la capacidad de resumen para comprobar los detalles del modelo:

summary(md_trend)

## 
## Call:
## lm(formula = y ~ trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5418 -0.3142 -0.1186  0.1732  2.7008 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.998e+00  3.673e-02 163.298   <2e-16 ***
## trend       7.359e-05  1.038e-04   0.709    0.479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4538 on 610 degrees of freedom
## Multiple R-squared:  0.0008228,  Adjusted R-squared:  -0.0008152 
## F-statistic: 0.5023 on 1 and 610 DF,  p-value: 0.4788

El valor proporcionado en la cuarta columna del último resultado de la regresión representa la probabilidad de que la hipótesis nula sea verdadera y rechazada, o de que el tipo I sea incorrecto. Los valores devueltos son 1, 0.1, 0.05, 0.01 y 0 Como siempre, proporciona contextualización de números a través de visualización de datos. Por lo tanto, usaremos el modelo diseñado para predecir el valor ajustado en la partición de entrenamiento y el valor predicho en la partición de prueba. La función de predicción del paquete de estadísticas, al igual que su nombre dado a entebder, predice el valor de los datos de entrada según el modelo dado. Lo usaremos para predecir los valores ajustados y predichos del modelo de tendencia que preparamos anteriormente.

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

Vamos a crear una función de utilidad que traza las series y la salida del modelo, utilizando el paquete plotly:

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

La premisa de la función es la siguiente:

Datos: Los datos de entrada, un objeto data.frame que sigue la misma estructura que la de Datos_sismo_df (incluyendo la variable yhat).
train: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo.
test: Igualmente, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de previsión.
title: El título de la parcela, por defecto, establecido en NULL.

Se constituye las entradas de la función plot_lm con la salida del modelo.

plot_lm(data = Datos_sismo_df,

train = train,

test = test,

title = "Predicción de Sismos para el componente de tendencia de la serie")

El modelo presentado fue capaz de captar el movimiento general de la tendencia, aunque es posible que una tendencia lineal no sea capaz de captar la ruptura estructural que se puede presentar en la tendencia alrededor de algunos de los años utilizados, como por ejemplo el caso evidenciado entre 1996 y 1998.

Para finalizar, para el análisis comparativo se busca medir la tasa de error del modelo tanto en el entrenamiento como en el conjunto de prueba, por eso se usa el código:

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

## [1] 0.05246603 0.07716665

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

Repasemos los detalles del modelo:

summary(md_seasonal)

## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6039 -0.3098 -0.1245  0.1686  2.5961 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.9725490  0.0637914  93.626   <2e-16 ***
## seasonalfeb   0.0003922  0.0902147   0.004    0.997    
## seasonalmar   0.0549020  0.0902147   0.609    0.543    
## seasonalabr   0.0960784  0.0902147   1.065    0.287    
## seasonalmay   0.1313725  0.0902147   1.456    0.146    
## seasonaljun   0.0960784  0.0902147   1.065    0.287    
## seasonaljul  -0.0333333  0.0902147  -0.369    0.712    
## seasonalago   0.0333333  0.0902147   0.369    0.712    
## seasonalsept  0.0490196  0.0902147   0.543    0.587    
## seasonaloct   0.0519608  0.0902147   0.576    0.565    
## seasonalnov   0.0588235  0.0902147   0.652    0.515    
## seasonaldic   0.0372549  0.0902147   0.413    0.680    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4556 on 600 degrees of freedom
## Multiple R-squared:  0.009441,   Adjusted R-squared:  -0.008719 
## F-statistic: 0.5199 on 11 and 600 DF,  p-value: 0.8905

Antes de trazar el modelo ajustado y los valores de previsión con la función plot_lm, vamos a actualizar los valores de yhat con la función predict, siendo el código:

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

Ahora podemos utilizar la función plot_lm para visualizar el modelo ajustado y los valores previstos, el código:

plot_lm(data = Datos_sismo_df,

train = train,

test = test,

title = "Prediccion de Sismos")

Como puede verse en el gráfico anterior, la tendencia de la serie no presenta una buen predicción. Antes de agregar los componentes de tendencia y estacional, evaluemos el desempeño 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.05247075 0.07604654

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

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

summary(md1)

## 
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5942 -0.3093 -0.1198  0.1663  2.6172 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.950e+00  7.112e-02  83.666   <2e-16 ***
## seasonalfeb   3.189e-04  9.025e-02   0.004    0.997    
## seasonalmar   5.476e-02  9.025e-02   0.607    0.544    
## seasonalabr   9.586e-02  9.025e-02   1.062    0.289    
## seasonalmay   1.311e-01  9.025e-02   1.452    0.147    
## seasonaljun   9.571e-02  9.025e-02   1.060    0.289    
## seasonaljul  -3.377e-02  9.025e-02  -0.374    0.708    
## seasonalago   3.282e-02  9.026e-02   0.364    0.716    
## seasonalsept  4.843e-02  9.026e-02   0.537    0.592    
## seasonaloct   5.130e-02  9.026e-02   0.568    0.570    
## seasonalnov   5.809e-02  9.026e-02   0.644    0.520    
## seasonaldic   3.645e-02  9.026e-02   0.404    0.686    
## trend         7.327e-05  1.043e-04   0.703    0.483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4558 on 599 degrees of freedom
## Multiple R-squared:  0.01026,    Adjusted R-squared:  -0.009571 
## F-statistic: 0.5173 on 12 and 599 DF,  p-value: 0.9042

El reingreso de la serie con la tendencia y los componentes estacionales juntos suministran un aumento adicional del R-cuadrado ajustado del modelo, que pasa de 0.004 a 0.02. Esta se puede observar en el gráfico de los resultados del modelo

train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = Datos_sismo_df,
train = train,
test = test,
title = "Prediccion de Sismos")

Midamos la puntuación MAPE del modelo en las particiones de entrenamiento y de prueba

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

## [1] 0.05236903 0.07633834

El uso de componentes de tendencia y estacionales para retroceder la serie puede mejorar significativamente la calidad y precisión del ajuste del modelo. Sin embargo, si observa el gráfico de ajuste y predicción del modelo, puede ver que la tendencia del modelo es demasiado lineal y no muestra un avance estructural en la tendencia de la serie. Aquí es donde agregar componentes polinomiales al modelo puede mejorar aún más la precisión del modelo. Una técnica sencilla para capturar tendencias no lineales es agregar un componente polinomial a la tendencia de secuencia para capturar la curvatura de la tendencia a lo largo del tiempo. Usamos el parámetro I para permitirnos realizar operaciones matemáticas en cualquier objeto de entrada. Entonces usamos este parámetro para agregar un polinomio cuadrático a la entrada de tendencia

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

El resumen del modelo puede verse como sigue:

summary(md2)

## 
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5959 -0.3077 -0.1204  0.1635  2.6236 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.943e+00  8.235e-02  72.163   <2e-16 ***
## seasonalfeb   3.177e-04  9.033e-02   0.004    0.997    
## seasonalmar   5.475e-02  9.033e-02   0.606    0.545    
## seasonalabr   9.586e-02  9.033e-02   1.061    0.289    
## seasonalmay   1.311e-01  9.033e-02   1.451    0.147    
## seasonaljun   9.571e-02  9.033e-02   1.060    0.290    
## seasonaljul  -3.378e-02  9.033e-02  -0.374    0.709    
## seasonalago   3.282e-02  9.033e-02   0.363    0.717    
## seasonalsept  4.843e-02  9.033e-02   0.536    0.592    
## seasonaloct   5.130e-02  9.033e-02   0.568    0.570    
## seasonalnov   5.809e-02  9.033e-02   0.643    0.520    
## seasonaldic   3.645e-02  9.033e-02   0.403    0.687    
## trend         1.480e-04  4.181e-04   0.354    0.724    
## I(trend^2)   -1.218e-07  6.605e-07  -0.184    0.854    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4561 on 598 degrees of freedom
## Multiple R-squared:  0.01031,    Adjusted R-squared:  -0.0112 
## F-statistic: 0.4793 on 13 and 598 DF,  p-value: 0.9363

train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = Datos_sismo_df,
train = train,
test = test,
title = "Predicción de tendencia (polinomial) para la clasificación de los Sismos")

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

## [1] 0.05235025 0.07620755

La función tslm

Hemos observado el proceso de transformación de un objeto ts a un formato de modelo de previsión de regresión lineal manualmente. La función tslm del paquete forecast ofrece una función integrada para transformar un objeto ts en un modelo de predicción de regresión lineal. Con la función tslm, puede establecer el componente de regresión junto con otras características.

Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie Datos_sismo con la función tslm utilizando la tendencia, el cuadrado de la tendencia y los componentes estacionales. En primer lugar, vamos a dividir la serie en particiones de entrenamiento y de prueba utilizando la función ts_split, cuyo código es:

Datos_sismo_split <- ts_split(Datos_sismo, sample.out = h)
train.ts <- Datos_sismo_split$train
test.ts <- Datos_sismo_split$test

A continuación, aplicaremos la misma fórmula que utilizamos para crear el modelo de previsión md2 anterior utilizando la función tslm:

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

Revisemos ahora md3, la salida de la función tslm, y comparémosla con la salida de md2:

summary(md3)

## 
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5959 -0.3077 -0.1204  0.1635  2.6236 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.943e+00  8.235e-02  72.163   <2e-16 ***
## season2      3.177e-04  9.033e-02   0.004    0.997    
## season3      5.475e-02  9.033e-02   0.606    0.545    
## season4      9.586e-02  9.033e-02   1.061    0.289    
## season5      1.311e-01  9.033e-02   1.451    0.147    
## season6      9.571e-02  9.033e-02   1.060    0.290    
## season7     -3.378e-02  9.033e-02  -0.374    0.709    
## season8      3.282e-02  9.033e-02   0.363    0.717    
## season9      4.843e-02  9.033e-02   0.536    0.592    
## season10     5.130e-02  9.033e-02   0.568    0.570    
## season11     5.809e-02  9.033e-02   0.643    0.520    
## season12     3.645e-02  9.033e-02   0.403    0.687    
## trend        1.480e-04  4.181e-04   0.354    0.724    
## I(trend^2)  -1.218e-07  6.605e-07  -0.184    0.854    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4561 on 598 degrees of freedom
## Multiple R-squared:  0.01031,    Adjusted R-squared:  -0.0112 
## F-statistic: 0.4793 on 13 and 598 DF,  p-value: 0.9363

Como puede observar en la salida anterior, ambos modelos (md2 y md3) son idénticos. El uso de la función tslm tiene varias ventajas, en lugar de establecer manualmente un modelo de regresión para las series con la función lm:

La eficiencia no requiere transformar las series a un objeto data.frame y ingeniería de características.
El objeto de salida soporta toda la funcionalidad de la previsión (como las funciones precisión y checkresiduals) y de los paquetes de TSstudio (como las funciones test_forecast y plot_forecast).

Modelización de eventos únicos y eventos no estacionales

En algunos casos, los datos de series de tiempo pueden contener o no patrones inusuales que se repiten con el tiempo. Los siguientes son ejemplos de estos tipos de eventos: - Valores atípicos: Un solo evento o eventos que están fuera del patrón normal de la serie. - Rotura estructural: hecho significativo que cambia los modelos históricos de la serie. Un ejemplo común es un cambio en el crecimiento de la serie. - Eventos no estacionales recurrentes: un evento que se repite de un ciclo a otro, pero su momento cambia de un ciclo a otro. Un ejemplo común de este tipo de eventos son las vacaciones de Semana Santa, que ocurren todos los años alrededor de marzo / abril.

Dado que este tipo de evento no se expresa en el modelo de regresión, sesgará los coeficientes estimados, ya que el modelo evaluará este tipo de evento junto con los eventos regulares de la serie.

El uso de variables de código activo, binarias o de marca podría ayudar al modelo a ignorar este tipo de eventos o ajustar los coeficientes del modelo en consecuencia.

En el diagrama de descomposición de la serie de gas de EE. UU. Que se muestra arriba, se puede ver, por ejemplo, que la tendencia de la serie alrededor de 2010 mostró una ruptura estructural. Si bien el crecimiento fue relativamente plano antes de 2010, la pendiente de la tendencia cambió posteriormente con un crecimiento positivo. En este caso, podemos usar una variable binaria que sea cero para las observaciones antes de 2010 y un año después.

La regresión de un modelo tslm con variables externas requiere un objeto data.frame separado con las variables correspondientes. El siguiente ejemplo muestra el proceso de creación de una variable binaria externa que es igual a 0 antes del año 2010 y a 1 después, utilizando tabla Datos_sismo_df

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

Ahora utilizaremos la nueva función para remodelar la serie Datos_sismo.

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

Utilicemos la función de resumen para revisar el resultado del modelo:

summary(md3)

## 
## Call:
## tslm(formula = Datos_sismo ~ season + trend + I(trend^2) + s_break, 
##     data = Datos_sismo_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6112 -0.3054 -0.1199  0.1773  2.6002 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.972e+00  8.457e-02  70.617   <2e-16 ***
## season2      2.594e-04  9.008e-02   0.003    0.998    
## season3      4.975e-02  9.008e-02   0.552    0.581    
## season4      7.078e-02  9.008e-02   0.786    0.432    
## season5      1.283e-01  9.008e-02   1.425    0.155    
## season6      8.398e-02  9.008e-02   0.932    0.352    
## season7     -6.037e-02  9.008e-02  -0.670    0.503    
## season8      1.834e-02  9.008e-02   0.204    0.839    
## season9      2.591e-02  9.008e-02   0.288    0.774    
## season10     2.097e-02  9.008e-02   0.233    0.816    
## season11     3.334e-02  9.009e-02   0.370    0.711    
## season12     1.590e-02  9.009e-02   0.177    0.860    
## trend       -3.731e-05  4.930e-04  -0.076    0.940    
## I(trend^2)   2.643e-07  8.771e-07   0.301    0.763    
## s_break     -1.257e-02  9.248e-02  -0.136    0.892    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4593 on 609 degrees of freedom
## Multiple R-squared:  0.01214,    Adjusted R-squared:  -0.01057 
## F-statistic: 0.5347 on 14 and 609 DF,  p-value: 0.9131

Como se puede observar en el resumen del modelo anterior, la variable de rotura estructural es estadísticamente significativa con un nivel de 0.01. La codificación en caliente también se puede aplicar a valores atípicos o días festivos estableciendo una variable binaria igual a 1 siempre que sea un valor atípico o un evento no estacional recurrente, 0.

De lo contrario, tenga en cuenta que después de haber entrenado un modelo de pronóstico con el La función tslm utiliza variables externas que deben generar valores futuros de estas variables, ya que se utilizarán como entrada para el pronóstico.

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

En los siguientes ejemplos usamos la serie UKgrid para demostrar la predicción de una serie de varias estaciones con un modelo de regresión lineal.

head(Datos_sismo_df)

##           ds   y trend seasonal s_break
## 1 1965-01-01 6.0     1      ene       0
## 2 1965-02-01 5.8     2      feb       0
## 3 1965-03-01 6.2     3      mar       0
## 4 1965-04-01 5.8     4      abr       0
## 5 1965-05-01 5.8     5      may       0
## 6 1965-06-01 6.7     6      jun       0

Como puede verse, esta serie tiene 4 variables: - ds: un objeto de fecha utilizado como timestamp o índice de series - 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) :

ts_heatmap (Datos_sismo ,  last  =  NULL ,  wday  =  TRUE ,  color  =  "greens" , 
  title  = "Mapeo de calor"  ,  padding  =  TRUE )

Podemos ver que es probable que observemos patrones de intensidad anual que abarcan 12 meses del año. Así, durante agosto de 1967, el evento más grande por magnitud se registró en la secuencia con una magnitud de 7,9. Los eventos de mayor intensidad sísmica se concentraron de marzo a octubre de los años 1965-1970, al igual que ocurrieron de 1993 a 2000 de muy baja intensidad sísmica.

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

Para obtener los componentes estacionales de la serie, la serie se añade 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)
Datos_sismo_dz<- Datos_sismo_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 Datos_sismo después de añadir estas nuevas características

str(Datos_sismo_dz)

## 'data.frame':    612 obs. of  8 variables:
##  $ ds      : Date, format: "1966-01-01" "1966-02-01" ...
##  $ y       : num  5.5 5.6 6 6.1 8.7 6 5.7 5.8 5.9 5.9 ...
##  $ trend   : int  13 14 15 16 17 18 19 20 21 22 ...
##  $ seasonal: Factor w/ 12 levels "ene","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 "dom\\."<"lun\\."<..: 7 3 3 6 1 4 6 2 5 7 ...
##  $ month   : Ord.factor w/ 12 levels "ene"<"feb"<"mar"<..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ lag365  : int  1 2 3 4 5 6 7 8 9 10 ...

Se usa la primera marca de tiempo de la serie y las funcionalidades year y yday (el día del año) del paquete lubridate para implantar el punto de inicio del objeto:

start_date <- min(Datos_sismo_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)
Datos_sismo_ts <- ts(Datos_sismo_dz$y, 
            start = c(year(start_date), yday(start_date)),
            frequency = 12
            )

Luego de cambiar la serie en un objeto ts, se puede volver atrás y confirmar la suposición que se logró sobre el grado de correlación entre la serie y sus intervalos estacionales con la funcionalidad ts_acf (una versión interactiva de la funcionalidad acf del paquete TSstudio). Se revisa la correlación de la serie con sus rezagos de los últimos cuatro años

library(TSstudio)
data(Datos_sismo_ts)

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

acf(Datos_sismo_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, luego de haber desarrollado las novedosas propiedades para el modelo y implantar el objeto ts, está todo listo para dividir la serie de acceso y el objeto de propiedades externas que corresponden que se inventaron (Datos_sismo) en una partición de entrenamiento y prueba. Como la finalidad es profetizar las próximas 365 visualizaciones, y la longitud de la serie es lo suficientemente enorme (2.540 observaciones), puede establecerse la partición de prueba a la longitud de las visualizaciones 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(Datos_sismo_ts, sample.out = h)
train_ts <- TRpartitions$train
test_ts <- TRpartitions$test

De igual manera, se separan las propiedades creadas para el modelo de regresión (las propiedades estacionales y de retardo) en una partición de entrenamiento y prueba siguiendo el mismo orden preciso utilizado para el objeto ts que corresponde.

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

Entrenamiento y prueba del modelo de previsión

Se usará la partición de entrenamiento y se empezarán a practicar los 3 modelos siguientes:

Modelo de alusión: Es la regresión de la serie con el elemento estacional y de tendencia en donde se usan las propiedades incorporadas de la funcionalidad tslm. Como se fija la frecuencia de la serie en 365, la característica estacional de la serie se refiere a la estacionalidad diaria.
Modelo multiestacional: En este modelo se agregan 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: Usando, además de los indicadores estacionales, la variable de retardo estacional.

Se comenzará con el modelo de alusión, en donde se va a hacer una regresión de la serie con sus elementos estacionales y de tendencia como se observa después:

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 = Datos_sismo_ts,
              forecast.obj = fc_tslm1,
              test = test_ts)

Se puede mirar en el gráfico de rendimiento anterior que el modelo de alusión hace un gran trabajo al capturar tanto la tendencia de la serie como la estacionalidad del día del año. Por otro lado, no consigue 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 funcionalidad 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 -2.223509e-17 0.4447976 0.3244619 -0.4884898 5.219197 0.7165313
## Test set      2.157769e-01 0.5871620 0.4934767  2.7108512 7.617690 1.0897780
##                    ACF1 Theil's U
## Training set 0.04993195        NA
## Test set     0.31748377 0.8031789

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(Datos_sismo)

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

Como puede verse en el gráfico de resumen de residuos anterior, los valores residuales no son blancos porque existe una correlación automática entre la serie residual y su rezago. Técnicamente, esto es una indicación de que el modelo no capturó todos los patrones o información en la secuencia. Una forma de abordar este problema es descubrir factores variables adicionales que permitan caracterizar la degradación de los sedimentos. El principal problema de este enfoque es la dificultad para detectar los factores externos predecibles que caracterizan la degradación de los sedimentos. Por ejemplo, los patrones sísmicos son difíciles de predecir con un mes de anticipación. Otra forma, cuando se dejan modelos en los residuos del modelo, es probarlos a partir de un modelo que consiste en probar los residuos del modelo como datos de series de tiempo independientes y modelarlos utilizando otro modelo de predicción de series de tiempo.

Resumen

Después de que cada proceso, comando, proceso y visualización se implementó en R, con datos de terremotos que ocurrieron con el trabajo en la tierra, fue posible analizar y organizar los conceptos clave desde el principio. Al inicio del informe dejar todo el resultado de las actividades y pedidos, claro sobre lo cumplido, solicitado y cumplido. Hemos entendido un poco más la importancia y organización de los datos geológicos y cómo interpretarlos correctamente en R.

Informe 3

Andres Garcia, Stefany Salazar, Angel Pacheco y Jhonatan Madrid

2021-11-04

PARTE 1

Características de la ingeniería de los componentes de la serie

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

La función tslm

Modelización de eventos únicos y eventos no estacionales

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

La serie UKgrid

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

Entrenamiento y prueba del modelo de previsión

Selección del modelo

Análisis de los residuos

Finalización de la previsión

Resumen

PARTE 2

Lectura de datos geologicos

Propiedades de la ingeniería de los componentes de la serie

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

La función tslm

Modelización de eventos únicos y eventos no estacionales

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

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

Entrenamiento y prueba del modelo de previsión

Análisis de los residuos

Resumen