Análisis y pronóstico de series de tiempo con R
Introducción
En este documento se explica lo relacionado al pronóstico de series de tiempo. Lo cual comprende:
Elaborando la serie de tiempo en R Análisis de la serie de tiempo Pronóstico de la serie de tiempo Estimación de errores de pronóstico Para el análisis de series de tiempo en R, existen una gran variedad de librerias. En este tutorial se emplea la libreria fpp2, la cual contiene contiene:
La libreria forecast, que contiene funciones de pronóstico La libreria ggplot2, que contiene funciones para gráficos Esta libreria se instala y se carga mediante install.packages(“fpp2”)y library(fpp2) respectivamente.
Elaborando la serie de tiempo en R Cargando librerias y datos Recordemos que se puede leer datos directamente desde una API, SQL u otros medios. Para este tutorial empleo la lectura de datos desde Excel. Entonces para ello emplearemos las siguientes funciones:
La función head() nos muestra los primeros datos de nuestra tabla de datos
library(fpp2)## Warning: package 'fpp2' was built under R version 4.2.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──## ✔ ggplot2 3.5.1 ✔ fma 2.5
## ✔ forecast 8.22.0 ✔ expsmooth 2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'forecast' was built under R version 4.2.3## Warning: package 'fma' was built under R version 4.2.3## Warning: package 'expsmooth' was built under R version 4.2.3## AirPassengers## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432Creando el objeto de serie de tiempo La serie de tiempo se debe almacenar en el objeto ts empleando la función ts(). El empleo de esta función es de la siguiente forma : serie <- ts(data, frequency= ,start=). Donde el primer agumento es un vector de datos. Luego los dos siguientes argumentos se indican de acuerdo a la siguiente tabla.
| Frecuencia | frequency= |
start= |
|---|---|---|
| Anual | 1 | 2000 |
| Trimestral | 4 | c(2000,2) |
| Mensual | 12 | c(2000,9) |
| Diario | 7 o 365.25 | 1 o c(2000,234) |
| Semanal | 52.18 | c(2000,23) |
| Horario | 24 o 168 o 8766 | 1 |
| Cada 30 min | 48 o 336 o 17532 | 1 |
En el siguiente fragmento de código definimos el objeto de la serie de tiempo y lo almacenaremos en data_serie.
#data_serie <- ts(data_serie$Serie, frequency=12, start=2004)
data_serie <- AirPassengers
data_serie## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432#Nota: Los resultados de esta guía corresponden a Data1 del archivo excel. Sugiero que tambien practiques con las otras series de tiempo como Data2 o Data3.
Análisis de la serie de tiempo Para el análisis de datos, esta libreria fpp2 nos permite emplear las siguientes funciones:
autoplot(), para graficar la serie de tiempo ggseasonplot(), para graficar la estacionalidad de una serie de tiempo ggsubseriesplot(), para graficar subseries ggAcf(), para graficar la autocorrelación decompose(), permite realizar una descomposición de estacionalidad y tendencia. Grafico de la serie de tiempo Para el gráfico simple de la serie de tiempo empleamos la funcion autoplot(). Donde solo se requiere ingresar el objeto de la serie de tiempo.
autoplot(data_serie)+
labs(title = "Serie de tiempo",
x = "Tiempo",
y = "Valor",
colour = "#00a0dc")+
theme_bw() 
Descomposición de la serie de tiempo Para la descomposición de la serie de tiempo se emplea la función decompose(). Donde se debe indicar el objeto de la serie de tiempo y el tipo de descomposición. Los tipos de descomposición que acepta esta función es additive y multiplicative.
Aditivo: Serie=T + S + I
Multiplicativo: Serie=T x S x I
Donde:
T: Tendencia
S: Estacionalidad
I: Irregular o error
# Descomposición de la serie de tiempo. Se almacena en el objeto fit
fit <- decompose(data_serie, type='additive')
#fit <- decompose(data_serie, type='multiplicative')
# Para graficar esta descomposición volvemos a emplear la funcion autoplot, pero con el objeto fit
autoplot(fit)+
labs(title = "Descomposicion de la serie de tiempo",
x = "Tiempo",
y = "Valor",
colour = "Gears")+
theme_bw()
library(highcharter)## Warning: package 'highcharter' was built under R version 4.2.3## Highcharts (www.highcharts.com) is a Highsoft software product which is## not free for commercial and Governmental usehchart(stl(data_serie, s.window='periodic'))## Warning: Deprecated function. Use the `create_axis` function.Grafico de la serie de tiempo con su tendencia El siguiente fragmento de código nos permite graficar la serie de tiempo con su tendencia. Notese que emplea el objeto fit en el cual guardamos previamente los valores de la descomposición. Nótese que se emplea la funcion trendcycle() para obtener los datos de tendencia del objeto fit.
autoplot(data_serie, series="Serie tiempo") +
autolayer(trendcycle(fit), series="Tendencia") +
labs(title = "Serie de tiempo",
x = "Tiempo",
y = "Valor"
) +
theme_bw()## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
Grafico de estacionalidad Para realizar el gráfico de estacionalidad empleamos la función ggseasonplot. Donde el argumento es el objeto que contiene la serie de tiempo.
ggseasonplot(data_serie)
library(TSstudio)## Warning: package 'TSstudio' was built under R version 4.2.3ts_seasonal(data_serie, type = "all")ts_heatmap(data_serie)#Graficar la autocorrelacion
acf(data_serie, lag=35, col="2",main ="",
xlab="Lag",ylab="Valores")
#Graficar la autocorrelacion parcial
pacf(data_serie, lag=35, col="2",main ="",
xlab="Lag",ylab="Valores")
#Como hay confirmación de tendencia se hace la diferencia (1 lag)
plot(diff(data_serie),xlab="Tiempo",ylab="Valores")
acf(diff(data_serie), lag=35, col="2",main = "",
xlab="Lag", ylab="Valores")
pacf(diff(data_serie), lag=35,col="2",main = "",
xlab="Lag",ylab="Valores")
#Como la variancia se ve incrementa en el tiempo se debe sacar logaritmo natural de la serie
plot(diff(log(data_serie)),xlab="Tiempo",ylab="Valores")
acf(diff(log(data_serie)), lag=35, col="2",main = "",
xlab="Lag", ylab="Valores")
pacf(diff(log(data_serie)), lag=35,col="2",main = "",
xlab="Lag",ylab="Valores")
#Augmented Dickey-Fuller test para ver no estacionaridad
#Devuelve la probabilidad de ser no estacionaria
#Si es mayor que 0.05 la serie es no estacionaria (tiene tendencia)
#Si es menor que 0.05 es estacionaria (No hay tendencia)
# library(aTSA)
aTSA::adf.test(diff(log(data_serie)))## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -9.61 0.01
## [2,] 1 -8.82 0.01
## [3,] 2 -7.63 0.01
## [4,] 3 -8.75 0.01
## [5,] 4 -6.79 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -9.63 0.01
## [2,] 1 -8.86 0.01
## [3,] 2 -7.71 0.01
## [4,] 3 -8.94 0.01
## [5,] 4 -6.98 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -9.60 0.01
## [2,] 1 -8.83 0.01
## [3,] 2 -7.69 0.01
## [4,] 3 -8.92 0.01
## [5,] 4 -6.95 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01library(forecast)
#Ajuste automático de los coeficientes
ARIMAmodel <- auto.arima(log(data_serie))
ARIMAmodel## Series: log(data_serie)
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4018 -0.5569
## s.e. 0.0896 0.0731
##
## sigma^2 = 0.001371: log likelihood = 244.7
## AIC=-483.4 AICc=-483.21 BIC=-474.77#Ajuste automático de los coeficientes
ARIMAmodel <- auto.arima(diff(log(data_serie)))
ARIMAmodel## Series: diff(log(data_serie))
## ARIMA(0,0,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4018 -0.5569
## s.e. 0.0896 0.0731
##
## sigma^2 = 0.001369: log likelihood = 244.7
## AIC=-483.39 AICc=-483.2 BIC=-474.77Pronóstico Métodos simples Para el pronóstico de series de tiempo mediante métodos básicos, la libreria fpp2 nos brinda las siguientes funciones:
naive(), metodo de naive simple ses(), exponential smoothing meanf(), media movil snaive(), metodo naive considerando estacionalidad El argumento a colocar en estas funciones es la serie de tiempo y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar. Por ejemplo si deseamos pronosticar 12 datos, se debe indicar h=12.
Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene un ajuste con la data historica checkresiduals(), permite analizar los residuales
# elaborando el método
m1 <- snaive(data_serie, h=96)
# graficando el pronóstico
autoplot(m1)
# verificando el ajuste del método
autoplot(m1)+autolayer(fitted(m1), series="Ajuste")## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
# verificando los residuales
checkresiduals(m1)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 275.04, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24Método regresión Para el pronóstico de series de tiempo mediante regresión, la libreria fpp2 nos brinda la función tslm(). Esta función la emplearemos para crear una regresión de la serie de tiempo con los datos de la descomposición estacional y/o tendencia. Entonces:
Si observamos solo tendencia usaremos : tslm(data_serie ~ trend) Si observamos solo estacionalidad usaremos : tslm(data_serie ~ season) Si observamos tendencia y estacionalidad usaremos : tslm(data_serie ~ trend + season) Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es la regresión y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.
Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene pronostico con la data historica checkresiduals(), permite analizar los residuales
library(tseries)## Warning: package 'tseries' was built under R version 4.2.3# elaborando la regresion
regresion <- tslm(data_serie ~ trend + season)
# elaborando el pronostico
m2 <- forecast(regresion)
# graficando el pronóstico
autoplot(m2)
# verificando el ajuste del método
autoplot(m2)+autolayer(fitted(m2), series="Ajuste")
# verificando los residuales
checkresiduals(m2)
##
## Ljung-Box test
##
## data: Residuals from Linear regression model
## Q* = 348.24, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24Método holt winters Para el pronóstico de series de tiempo mediante holt winters, la libreria fpp2 nos brinda la función hw().
Los argumentos a colocar en esta funcion son:
La serie de tiempo El valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar. El tipo de descomposición a usar para la estacionalidad. Los tipos de descomposición que acepta esta función es additive y multiplicative. Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene pronostico con la data historica checkresiduals(), permite analizar los residuales
# elaborando el pronostico
m3 <- hw(data_serie, h=96, seasonal = 'multiplicative')
# graficando el pronóstico
autoplot(m3)
# verificando el ajuste del método
autoplot(m3)+autolayer(fitted(m3), series="Ajuste")
# verificando los residuales
checkresiduals(m3)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 40.151, df = 24, p-value = 0.0206
##
## Model df: 0. Total lags used: 24ARIMA Para el pronóstico de series de tiempo mediante ARIMA, la libreria fpp2 nos brinda la función auto.arima().
Primero crearemos un modelo ARIMA, para ello el argumento a colocar en esta funcion es la serie de tiempo. Considerar que esta función es solo una aproximación iterativa que busca los indices de AR y MA. Pues en determinadas series de tiempo podria no encontrar los indices adecuados para un modelo ARIMA. En ese caso lo adecuado es seguir la metodología de estimación de índices ARIMA. Esta metodología no esta cubierta en esta guia.
Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es el modelo ARIMA y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.
Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene pronostico con la data historica checkresiduals(), permite analizar los residuales
# elaborando el modelo ARIMA
modelo_arima <- auto.arima(data_serie)
# elaborando el pronostico
m4 <- forecast(modelo_arima, h=96)
# graficando el pronóstico
autoplot(m4)
# verificando el ajuste del método
autoplot(m4)+autolayer(fitted(m4), series="Ajuste")
# verificando los residuales
checkresiduals(m4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)(0,1,0)[12]
## Q* = 37.784, df = 21, p-value = 0.01366
##
## Model df: 3. Total lags used: 24Red neuronal Para el pronóstico de series de tiempo mediante una red neuronal, la libreria fpp2 nos brinda la función nnetar().
Primero crearemos un modelo de red neuronal (neural network), para ello el argumento a colocar en esta funcion es la serie de tiempo.
Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es el modelo de red neuronal y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.
Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene pronostico con la data historica checkresiduals(), permite analizar los residuales
# elaborando el modelo de red neuronal
neural_network <- nnetar(data_serie)
# elaborando el pronostico
m5 <- forecast(neural_network, h=96)
# graficando el pronóstico
autoplot(m5)
# verificando el ajuste del método
autoplot(m5)+autolayer(fitted(m5), series="Ajuste")## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
# verificando los residuales
checkresiduals(m5)
##
## Ljung-Box test
##
## data: Residuals from NNAR(1,1,2)[12]
## Q* = 195.94, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24Estimación de error Para estimar los errores de pronóstico, se debe realizar con los valores ocurridos o reales. Es decir este error se mide tiempo despues de haber realizado el pronóstico.
Supongamos que los valores reales ocurridos para la Data1 son los siguientes:
Mes Valor Real
Enero 19 13487 Mayo Febrero 19 12776 Junio Marzo 19 13812 Julio Abril 19
13032 Agosto Entonces, para poder comparar nuestros datos sera necesario
almacenarlo en un objeto de serie de tiempo. Entonces:
real <- c(622, 606, 508, 461, 390, 432,462, 471)
data_real <- ts(real, frequency=12,start=1961)La función accuracy() determina los errores de pronostico, para ello es necesario:
Indicar el modelo de pronostico, donde estará el pronóstico Indicar datos reales, donde estará los valores reales Entonces lo que realizará esta función es comparar el pronóstico de los siguientes 8 datos y el valor real. Pues estamos considerando que ya pasaron 8 meses, y nos encontramos en la etapa de evaluar el error de pronóstico de nuestros modelos.
# modelo en base a métodos simples
accuracy(m1,data_real)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 31.77273 36.31574 32.0303 11.124393 11.24871 1.000000 0.7464603
## Test set 3.62500 140.23596 123.6250 -2.775782 24.26489 3.859626 0.6988219
## Theil's U
## Training set NA
## Test set 2.475277# modelo en base a regresion lineal
accuracy(m2,data_real)## ME RMSE MAE MPE MAPE MASE
## Training set 6.123044e-15 25.11363 19.77364 0.6039558 8.591515 0.6173418
## Test set -9.952652e-01 106.24590 92.63684 -3.1086876 18.289105 2.8921624
## ACF1 Theil's U
## Training set 0.7631551 NA
## Test set 0.6491973 1.800463# modelo en base a holt winters
accuracy(m3,data_real)## ME RMSE MAE MPE MAPE MASE
## Training set 1.256973 10.63256 7.790649 0.2182707 2.914411 0.2432275
## Test set -35.415680 150.14883 136.598759 -10.9560051 27.796830 4.2646727
## ACF1 Theil's U
## Training set 0.2135914 NA
## Test set 0.6873988 3.013635# modelo en base a ARIMA
accuracy(m4,data_real)## ME RMSE MAE MPE MAPE MASE
## Training set 1.34230 10.84619 7.86754 0.420698 2.800458 0.245628
## Test set -27.27284 144.07536 132.46579 -9.203000 26.842386 4.135640
## ACF1 Theil's U
## Training set -0.00124847 NA
## Test set 0.69852717 2.839455# modelo en base a red neuronal
accuracy(m5,data_real)## ME RMSE MAE MPE MAPE MASE
## Training set 0.006053123 14.47562 11.59752 -0.2947613 4.273755 0.3620796
## Test set -47.260966160 157.21117 144.00954 -13.4509792 29.632264 4.4960406
## ACF1 Theil's U
## Training set 0.5604998 NA
## Test set 0.7016088 3.274349Entonces a partir de estos errores de pronóstico podemos determinar cual es el modelo mas adecuado. En la práctica se suelen evaluar varios modelos e incluso tomar como pronóstico un valor medio del resultados de dos o más modelos.