INTRODUCCIÓN

El presente informe, una vez más, nos convence de la importancia de el análisis estadístico en la cotidianidad y en la resolución de problemáticas actuales y reales, que se desarrollan a nuestro alrededor y que nos involucran. Este trabajo, pretende mostrar cómo a través del servidor de R analizar, identificar y practicar cada uno de los siguientes temas: 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, apreciació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 UK daily, Entrenamiento y prueba del modelo de previsión, selección del modelo, Analisis de residuos, Finalización de previsión.

Parte 1

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.Primero debemos entender su composición. En los siguientes ejemplos se muestra el proceso de construcción de nuevas propiedades a partir de lotes existentes utilizando usando la serie USgas.

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

El resultado se muestra a continuacion:

library(TSstudio)

data(USgas)

ts_plot(USgas,

title =
"US Monthly Natural Gas consumption",

Ytitle =
"Billion Cubic Feet",

Xtitle =
"Year")

Igualmente,revisaremos características primordiales de la serie utilizando la función ts_info, donde su resultado es el siguente:

ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

En los capitulos anteriores USgas 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,como se muestra a continuacion el resultado:

ts_decompose(USgas)

La tendencia de la serie es bastante plana entre 2000 y 2010 en el grafico mostrado anteriormente ,y tiene un aumento bastante lineal hacia adelante. Por eso la tendencia general entre 2000 y 2018 no es exactamente lineal. Este es un dato esencial que nos permitirá conceptualizar el ingreso de la tendencia para el modelo de regresión.

Se tiene que cambiar la serie de un objeto ts a un objeto data.frame para modificar la capacidad de regresión lineal incorporada en R del paquete stats antes de usar la funcionalidad lm,

Por eso,vamos a usar la funcionalidad ts_to_prophet del paquete TSstudio, que tiene como resultado:

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. El código es el siguiente:

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

Despues de modificar las series en un objeto data.frame, tenemos la posibilidad de empezar a generar las propiedades de acceso de la regresión. La primera característica que crearemos es la tendencia de la serie. Un enfoque que no se puede sustituir o cambiar para edificar la variable de tendencia y hacer indices de las visualizaciones de la serie en orden cronológico, el código reultante seria:

USgas_df$trend <- 1:nrow(USgas_df)

La regresión de la serie con el índice de la serie da una valoracion del incremento secundario mensual, debido a que el índice está en orden cronológico con incrementos constantes.

La segunda característica que se desea elaborar es el elemento estacional. Ya que deseamos medir el aporte 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 simbolizan los meses del año y, por consiguiente, crearemos una variable categórica con 12 categorías,las cuales corresponde a un mes del año.Ademas, Usaremos la funcionalidad month del paquete lubridate para sustraer el mes del año de la variable de fecha ds, el código coligado es:

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)

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, el resultado para el código es:

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

Previo a comenzar a hacer la regresión de la serie con dichas propiedades, se 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), ]

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

El coeficiente de la variable de tendencia es estadísticamente primordial a un grado de 0,001.Sin embargo, 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.

La cuarta columna del resultado de la regresion anterior representa el costo p de todos los coeficientes del modelo. El costo p concede la posibilidad de que se rechace la hipótesis nula ya que es verdadera, o el error de tipo I. Entonces, para el costo p menor que, el costo umbral, rechazaremos la hipótesis nula con un grado de importancia, donde los valores normales de son 0,1, 0,05, 0,01, etcétera.

Como constantemente, se ofrece contextualizar los números con la visualización de los datos. Por lo tanto, usaremos el modelo que se ha diseñado para predecir 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 da a enteder, predice los valores de un dato de ingreso basándose en un modelo ya dado.

La usaremos para predecir los valores ajustados y previstos del modelo de tendencia que hemos preparado antes, el código es:

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)

## 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)
}

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

Se constituye las entradas de la función plot_lm con la salida del modelo. El código es:

plot_lm(data = USgas_df,

train = train,

test = test,

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

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

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.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, el código es el siguiente:

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, 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 = 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, código es:

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:l, el código asociado es:

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 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,78 a 0,91. Esta se puede observar en el gráfico de los resultados del modelo, el código es:

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, siendo el código:

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, cuyo código asociado es:

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%, el código es:

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

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 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, cuyo código es:

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)

## 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:

*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 USgas_df, cuyo código es:

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. El código asociado es:

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

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

summary(md3)

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

Como 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. El código asociado es:

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

Utilizaremos la función head para revisar las variables de la serie, cuyo código es:

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, el código es:

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

-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, el código asociado es:

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, cuyo código asociado es:

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 el programa tiene una relación sólida durante el fin de semana (es decir, 7, 1, 21, etc.). Por lo tanto, en la mayoría de los casos, no se recomienda usar minutos cortos en lugar de predecir (por ejemplo, la séptima vez o la hipótesis es 365) porque debes usar estos momentos como un regalo en la profecía de la naturaleza de la profecía. Esto requiere un cierre adicional porque tiene que predecir estos retrasos. Además, puede agregarnos variaciones predictivas utilizando métodos predictivos como la ubicación.

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

0 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, el código es:

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 GRT. 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, el código es:

md_tslm1 <- tslm(train_ts ~ season + trend)

Ahora,utilizaremos el modelo entrenado md_tslm1 para pronosticar los próximos 365 días seguidos teniendo en cuenta las observaciones de la partición de prueba con el predictor:

fc_tslm1 <- forecast(md_tslm1, h = h)

Al comparar 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)

Podemos ver que el índice de referencia hace un buen trabajo al apresar 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.Desde otro punto de vista, no se registra la vibración relacionada con el día de la semana.Se puede ver en el resultado de la siguiente función de precisión, 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

Se intenta mejorar la precisión del modelo, sumandole al modelo las características del día de la semana y del mes del año. El código es:

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

Como ahora estamos utilizando características de una fuente de datos externa, tenemos que especificar 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, cuyo código asociado:

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)

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

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, el código asociado es:

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 resumen para observar el nivel de significancia de la variable lag3655. El código es: Codigo aquí:

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.

Backtest: El tercer modelo puede ser correcto solo por casualidad, no porque la varianza contribuya a la precisión del modelo, ya que la diferencia es 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.

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, los modos se analizan con equipo de entrenamiento, incluido el procedimiento seleccionado. Para simplificar, nos atenemos a los criterios de precisión y elegimos el tercer modelo para predecir la serie. El siguiente paso es reconstruir el modelo para todas las filas y predecir 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 puede verse en el resumen anterior de residuos, los residuos no son ruido blanco porque existe cierta autocorrelación entre la serie residual y su intervalo. Esta es una indicación técnica de que el modelo no contiene todos los patrones o la información de la serie. Una forma de resolver este problema es identificar variables adicionales que puedan explicar las variaciones en los residuos. El principal problema con este enfoque es que es difícil identificar variables atípicas que puedan explicar y predecir las variaciones residuales. 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 objeto data.frame con los valores de wday, month y lag365 para las próximas 365 observaciones futuras. Un enfoque simplificado es establecer primero las observaciones pronosticadas y luego extraer el día de la semana y el mes del año de este objeto. El código asociado es:

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

Podemos utilizar la variable date para establecer las variables wday y month con el paquete lubridate como se muestra a continuacion:

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 de previsión para crear la previsión. El código coligado es:

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 presentamos las aplicaciones predictivas de los modelos de regresión lineal. La principal ventaja de los modelos de regresión lineal sobre otros modelos tradicionales de series de tiempo es la capacidad del modelo para tener en cuenta variables y valores atípicos. Sin embargo, este modelo puede manejar series de tiempo con el modelo de estaciones múltiples, como vimos con el pronóstico de demanda de electricidad en el Reino Unido..

Parte 2

Primeramente,debemos cargar los datos del archivo CSV a R haciendo una ruta para ubicar el archivo que se quiere utilizar, la ruta consiste en dos partes: *Primera parte: nombre del archivo. se debe asegurar que el nombre sea idéntico al nombre del archivo real que se va a cargar *segunda parte: tipo de archivo.csv

library(readr)
Datos_geologicos <- read.csv("C:/Users/albon/Desktop/database.csv")
class(Datos_geologicos)

## [1] "data.frame"

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

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

class(Datos_geologicos)

## [1] "ts"

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

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

El resultado se muestra a continuacion:

library(TSstudio)

data(Datos_geologicos)

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

ts_plot(Datos_geologicos,

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, donde su resultado es el siguente:

ts_info(Datos_geologicos)

##  The Datos_geologicos 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_geologicos 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,como se muestra a continuacion el resultado:

ts_decompose(Datos_geologicos)

La tendencia de la serie estacional es creciente a medida que avanzan los meses del año especificado hasta llegar a su maximo valor de magnitud, a lo cual empezara a descender hasta empezar el siguente año, esto se aplica en los 53 años presentados.Por eso, la tendencia general entre 1965 y 2016 no es lineal. Este es un dato esencial que nos permitirá conceptualizar el ingreso de la tendencia para el modelo de regresión.

Se tiene que cambiar la serie de un objeto ts a un objeto data.frame para modificar la capacidad de regresión lineal incorporada en R del paquete stats antes de usar la funcionalidad lm,

Por eso,vamos a usar la funcionalidad ts_to_prophet del paquete TSstudio, que tiene como resultado:

Datos_geologicos_df <- ts_to_prophet(Datos_geologicos)

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

Datos_geologicos_df$trend <- 1:nrow(Datos_geologicos_df)

La regresión de la serie con el índice de la serie da una valoracion del incremento secundario mensual, debido a que el índice está en orden cronológico con incrementos constantes.

La segunda característica que se desea elaborar es el elemento estacional. Ya que deseamos medir el aporte 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 Datos_geologicos, las unidades de frecuencia simbolizan los meses del año y, por consiguiente, crearemos una variable categórica con 12 categorías,las cuales corresponde a un mes del año.Ademas, Usaremos la funcionalidad month del paquete lubridate para sustraer el mes del año de la variable de fecha ds, el código coligado es:

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

head(Datos_geologicos_df)

##           ds   y trend seasonal
## 1 1965-01-01 6.5     1      ene
## 2 1965-02-01 5.9     2      feb
## 3 1965-03-01 5.6     3      mar
## 4 1965-04-01 5.6     4      abr
## 5 1965-05-01 6.0     5      may
## 6 1965-06-01 6.8     6      jun

h <- 12 # setting a testing partition length
train <- Datos_geologicos_df[1:(nrow(Datos_geologicos_df) - h), ]
test <- Datos_geologicos_df[(nrow(Datos_geologicos_df) - h + 1):nrow(Datos_geologicos_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.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

El valor proporcionado en la cuarta columna del resultado de la regresion anterior representa la posibilidad de que se rechace la hipótesis nula ya que es verdadera, o el error de tipo I.Donde los valores arrojados son son 1, 0,1, 0,05, 0,01 y 0

Como constantemente, se ofrece contextualizar los números con la visualización de los datos. Por lo tanto, usaremos el modelo que se ha diseñado para predecir 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 da a entebder, predice los valores de un dato de ingreso basándose en un modelo ya dado.

La usaremos para predecir los valores ajustados y previstos del modelo de tendencia que hemos preparado antes, el código es:

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_geologicos_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. El código es:

plot_lm(data = Datos_geologicos_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.05435221 0.04106426

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.4980 -0.3054 -0.1569  0.1657  1.9588 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.87255    0.06184  94.959   <2e-16 ***
## seasonalfeb  -0.01569    0.08746  -0.179    0.858    
## seasonalmar   0.09216    0.08746   1.054    0.292    
## seasonalabr   0.03922    0.08746   0.448    0.654    
## seasonalmay   0.08902    0.08746   1.018    0.309    
## seasonaljun   0.12549    0.08746   1.435    0.152    
## seasonaljul   0.03725    0.08746   0.426    0.670    
## seasonalago   0.06863    0.08746   0.785    0.433    
## seasonalsept -0.05882    0.08746  -0.673    0.501    
## seasonaloct   0.11569    0.08746   1.323    0.186    
## seasonalnov  -0.06863    0.08746  -0.785    0.433    
## seasonaldic  -0.06078    0.08746  -0.695    0.487    
## ---
## 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

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_geologicos_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, código es:

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

## [1] 0.05383805 0.04855883

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

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.5457 -0.2986 -0.1457  0.1629  1.9132 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.9223487  0.0688236  86.051   <2e-16 ***
## seasonalfeb  -0.0155208  0.0873361  -0.178    0.859    
## seasonalmar   0.0924878  0.0873363   1.059    0.290    
## seasonalabr   0.0397120  0.0873366   0.455    0.649    
## seasonalmay   0.0896814  0.0873370   1.027    0.305    
## seasonaljun   0.1263174  0.0873375   1.446    0.149    
## seasonaljul   0.0382476  0.0873382   0.438    0.662    
## seasonalago   0.0697856  0.0873389   0.799    0.425    
## seasonalsept -0.0574999  0.0873398  -0.658    0.511    
## seasonaloct   0.1171753  0.0873408   1.342    0.180    
## seasonalnov  -0.0669730  0.0873419  -0.767    0.444    
## seasonaldic  -0.0589644  0.0873431  -0.675    0.500    
## trend        -0.0001654  0.0001009  -1.639    0.102    
## ---
## 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

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, el código es:

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

Midamos la puntuación MAPE del modelo en las particiones de entrenamiento y de prueba, siendo el código:

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)

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.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 ***
## seasonalfeb  -1.550e-02  8.683e-02  -0.179  0.85836    
## seasonalmar   9.252e-02  8.683e-02   1.066  0.28707    
## seasonalabr   3.976e-02  8.683e-02   0.458  0.64723    
## seasonalmay   8.973e-02  8.683e-02   1.033  0.30184    
## seasonaljun   1.264e-01  8.683e-02   1.455  0.14610    
## seasonaljul   3.830e-02  8.683e-02   0.441  0.65930    
## seasonalago   6.984e-02  8.683e-02   0.804  0.42157    
## seasonalsept -5.746e-02  8.683e-02  -0.662  0.50843    
## seasonaloct   1.172e-01  8.684e-02   1.350  0.17760    
## seasonalnov  -6.696e-02  8.684e-02  -0.771  0.44098    
## seasonaldic  -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

train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = Datos_geologicos_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.05325876 0.05659700

La función tslm

Ahora repetiremos el ejemplo anterior y pronosticaremos las últimas 12 observaciones de la serie Datos_geologicos 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_geologicos_split <- ts_split(Datos_geologicos, sample.out = h)
train.ts <- Datos_geologicos_split$train
test.ts <- Datos_geologicos_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.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 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.

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.

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

Ahora utilizaremos la nueva función para remodelar la serie Datos_geologicos. El código asociado es:

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

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

summary(md3)

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

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

##           ds   y trend seasonal s_break
## 1 1965-01-01 6.5     1      ene       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      abr       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) :

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

Se puede observar que tenemos la posibilidad de mirar patrones de intensidad anuales distribuidos en los 12 meses del año .Es asi como para el mes de Agosto del año 1967 se registra el más grande acontecimiento de intensidad de la serie con intensidad de 7,9. 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 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_geologicos_dz<- Datos_geologicos_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_geologicos después de añadir estas nuevas características, obteniendo la siguiente salida:

str(Datos_geologicos_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: 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_geologicos_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_geologicos_ts <- ts(Datos_geologicos_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, obteniendo la gráfica de salida:

library(TSstudio)
data(Datos_geologicos_ts)

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

acf(Datos_geologicos_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_geologicos) 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_geologicos_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_geologicos_df[1:(nrow(Datos_geologicos_df) - h), ] 
test_df <- Datos_geologicos_df[(nrow(Datos_geologicos_df) - h + 1):nrow(Datos_geologicos_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_geologicos_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  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(Datos_geologicos)

## 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 sonido blanco, ya existente cierta autocorrelación en medio de las series de residuos y sus rezagos. Es decir 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 detectar cambiantes extras que logren describir la alteración de los residuos. El primordial problema de este enfoque es que es complicado detectar cambiantes externas que logren describir la alteración de los residuos y que sean factibles de adivinar. Ejemplificando, es complicado presagiar patrones sismicos con un mes de antelación. Un enfoque alternativo, una vez que los patrones sobran en los residuos del modelo, es intentarlos de un modelo es intentar los residuos del modelo como datos de una serie temporal sin dependencia y modelizarlos con otro modelo de previsión de series temporales.

CONCLUSIÓN

Luego de realizar cada uno de los procesos, mandos, operaciones y visualizaciones en R, de los temas mencionados en la introducción de este informe, y con datos de sismos ocurridos con los cuales se trabajo una Marte del trabajo, se pudo analizar y organizar los conceptos base al inicio del informe, dejando cada uno de los resultados de operaciones y mandos, claridad acerca de lo ejecutado, ordenado y finalizado. Comprendimos un poco más la importancia y organización de datos geológicos y como interpretarlos en el servidor de R de manera correcta, usando nuevos caracteres y algoritmos para ejecutar las funciones requeridas y solicitadas desde principios de entrega del informe de Análisis Estadístico

Informe #3

Armando Bonilla, Camila Simanca, Jaiza Causil, Camila Vasquez

1/11/2021

INTRODUCCIÓN

Parte 1

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

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

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

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

CONCLUSIÓN