Introducción

Una serie de tiempo es una colección de observaciones de una variable tomadas de forma secuencial y ordenada en el tiempo (instantes de tiempo equiespacios). Las series pueden tener periodicidad anual, semestral, trimestral, mensual, etc., según los periodos de tiempo en los que están recogidos los datos que la componen. Las series temporales se pueden definir como un caso particular de los procesos estocásticos, ya que un proceso estocástico es una secuencia de variables aleatorias, ordenadas y equidistantes cronológicamente referidas a una características observables en diferentes momentos.

Algunos ejemplos de series temporales vienen de campos como la economía (producto interior bruto anual, tasa de inflación, tasa de desempleo, …), la demografía (nacimientos anuales, tasa de dependencia, …), la meteorología (temperaturas máximas, medias o mínimas, precipitaciones diarias, …), etc.

El análisis de series temporales explica el hecho de que los puntos de datos tomados a lo largo del tiempo pueden tener una estructura interna (como la autocorrelación, la tendencia o la variación estacional) que debe tenerse en cuenta.http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc41.htm

Aplicaciones

Obtener una comprensión de las fuerzas subyacentes y la estructura que produjo los datos observados
Ajustar un modelo y proceder a la previsión, el control o incluso la retroalimentación y el control anticipativo.

El análisis de series de tiempo se utiliza para muchas aplicaciones tales como:

Previsión económica
Pronóstico de ventas
Análisis presupuestario
Análisis del mercado de valores
Proyecciones de rendimiento
Proceso y control de calidad
Estudios de inventario
Proyecciones de carga de trabajo
Estudios de utilidad
Análisis del censo

Ejemplo

Objetivo

El objetivo es el análisis, modelado y predicción de la serie temporal CO2

Análisis de la Serie

Utilizaremos co2 que puedes descargar del siguiente enlace https://www.dropbox.com/s/zugg2443hlq7z6z/co2.csv?dl=0

#Se carga las librerias necesarias
library(forecast)
library(tseries)

#Cargamos los datos co2
co2<-read.csv("C:/datos/co2.csv", header = TRUE, sep = ";")
#Transformamos los datos en una serie temporal 
co2ts<-ts(co2$CO2, start = c(1959,1), frequency = 12)
print(co2ts)

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
## 1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
## 1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
## 1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
## 1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
## 1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
## 1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
## 1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
## 1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
## 1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
## 1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
## 1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
## 1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
## 1972 351.71 350.94 349.10 346.77 345.73                                   
##         Nov    Dec
## 1959 331.87 333.18
## 1960 333.57 334.72
## 1961 334.86 336.07
## 1962 336.54 337.79
## 1963 337.93 338.95
## 1964 340.10 340.88
## 1965 341.69 342.51
## 1966 342.70 343.65
## 1967 343.04 345.27
## 1968 345.11 347.07
## 1969 347.23 347.62
## 1970 347.38 349.38
## 1971 349.29 350.91
## 1972

library(ggfortify)

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.4.2

## 
## Attaching package: 'ggfortify'

## The following object is masked from 'package:forecast':
## 
##     gglagplot

#Trazamos la serie de tiempo co2ts
autoplot(co2ts, ts.colour = "blue", ts.linetype = "dashed")

Trazada ya la serie co2ts se puede determinar que presenta una tendencia ascendente (la serie no es estacionaria en media). Con esto se puede afirmar que la serie no es estacionaria. Tambien se puede observar que la serie es estacionaria en cuando a la varianza, ya que no se aprecia gran variabiliad. Para confirmar lo dicho vamos a trazar la gráfica de su función de autocorrelación.

autoplot(acf(co2ts, plot = FALSE))

Se puede apreciar en esta gráfica como la serie no es estacionaria ya que el valor de la función de autocorrelación no decae de manera exponencial a medida que aumentan los rezagos en el tiempo. También apreciamos que existe la componente estacional, con un periodo de 12, por lo tanto se analizará las componente por separado para ello descomponemos la serie.

autoplot(stl(co2ts, s.window = "periodic"), ts.colour = "blue")

Modelado de la Serie

El análisis previo nos revela que tenemos que eliminar la tendencia y la estacionalidad de la serie para que pueda ser estacionaria. El modelo que vamos a utilizar es un modelo ARIMA, diferenciando la serie lograremos que sea convertida a estacionaria, por lo tanto vamos a ver que nos dice las funciones ndiffs y nsdiffs, que calculan cada una el número de diferenciaciones regulares y estacionales respectivamente que se necesita llevar a cabo para que la serie sea estacionaria.

ndiffs(co2ts)

## [1] 1

nsdiffs(co2ts)

## [1] 1

Los respectivos cálculos nos muestra que la serie necesitar ser diferenciación regular y otra estacional.

diff.co2ts<-autoplot(diff(co2ts), ts.linetype = "dashed", ts.colour = "darkmagenta")
diff.co2ts

Como podemos ver ahora en la serie ya se ha eliminado la componente de tendencia, una ves eliminada esta componente volvemos a gráficar la función de autocorrelación.

autoplot(acf(diff(co2ts), plot = FALSE))

Podemos apreciar que la serie sigue sin ser estacionaria, ya que la función de autocorrelación decrece rápidamente en los desfaces regulares de forma lenta en los retardos (12,24,36), de forma que no es estacionaria en la parte estacional, además nuevamente se vuelve a ver la componente estacional de perido 12.

monthplot(diff(co2ts), col = "midnightblue")

diffco2<-diff(co2ts)
boxplot(diffco2~cycle(diffco2))

Se puede apreciar que los picos más altos se presentan en los meses de marzo y abril, mientras que los picos más bajos se presentan en los meses agosto y septiembre. Ahora eliminaremos el componente estacional

diff.co2ts.12<-diff(co2ts, lag = 12)
autoplot(diff.co2ts.12, ts.colour = "darkorange4", ts.linetype = "dashed")

Una vez eliminado tanto la componente de tendencia y la estacional, observemos que esta serie se parece bastante a una serie estacionaria, ya que parece ser constante en media y varianza, para asegurarnos aplicamos un test de estacionariedad como el test ADF (Dickey-Fuller) y el test de KPSS (Kwiatkowski-Phillips-Schmidt-Shin).

Planteamiento de Hipótesis
Significancia \(\alpha\) = 0.05

\(H_{0}:\) La serie es no estacionaria: tiene raíz unitaria

\(H_{1}:\) La serie es estacionaria: no tiene raíz unitaria

library(tseries)
adf<-adf.test(diff.co2ts.12)
adf$p.value

## [1] 0.02201398

kpss<-kpss.test(diff.co2ts.12)

## Warning in kpss.test(diff.co2ts.12): p-value greater than printed p-value

kpss$p.value

## [1] 0.1

El valor 0.01 del test indica que se puede rechazar la hipótesis nula \(H_{0}\) indicando que la serie es es no estacionaria.

autoplot(acf(diff.co2ts.12, plot = FALSE))

autoplot(pacf(diff.co2ts.12, plot = FALSE))

Vistas la gráficas acf y pacf pueden plantear varios modelos tentativos para el análisis, como lo serían los ARIMA(0,1,2)(0,1,1) o ARIMA(1,1,0)(2,1,0) deducidos, o realizando una combinación de estos como ARIMA(1,1,2)(2,1,1)) o incluso variando algunos de los componentes como puede ser con ARIMA(1,1,1)(2,1,1)) o ARIMA(1,1,2)(1,1,1)) o ARIMA(0,1,1)(0,1,1) o ARIMA(1,1,0)(1,1,0), etc.

library(forecast)
arima1<- Arima(co2ts, order=c(0,1,2), seasonal=list(order=c(0,1,1),period=12))
arima2<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(2,1,0),period=12))
arima3<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(2,1,1),period=12))
arima4<- Arima(co2ts, order=c(1,1,1), seasonal=list(order=c(2,1,1),period=12))
arima5<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(1,1,1),period=12))
arima6<- Arima(co2ts, order=c(0,1,1), seasonal=list(order=c(0,1,1),period=12))
arima7<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(1,1,0),period=12))

Cálculamos los criterios AIC y BIC.

AIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)

##        df       AIC
## arima1  4  85.72537
## arima2  4  93.70131
## arima3  7  90.55929
## arima4  6  89.02326
## arima5  6  88.62985
## arima6  3  83.82626
## arima7  3 106.19086

BIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)

##        df       BIC
## arima1  4  97.71422
## arima2  4 105.69016
## arima3  7 111.53978
## arima4  6 107.00654
## arima5  6 106.61312
## arima6  3  92.81789
## arima7  3 115.18250

Aquí se puede apreciar que los ajustes que mejor AIC y BIC presentan son aquellas que solo tienen componente de médias móviles y no tienen componente autorregresiva, siendo ARIMA(0,1,1)(0,1,1) el modelo que los test arrojan con un menor valor y por tanto con una mayor consideración. Una vez estimados los modelos y elegido el mejor de ellos, en este caso ARIMA(0,1,1)(0,1,1), se procede a validarlo. Para ello en primer lugar se grafican los correlogramas de los residuos para comprobar que son ruido blanco:

autoplot(acf(arima6$residuals, plot = FALSE))

autoplot(pacf(arima6$residuals, plot = FALSE))

Se puede apreciar en los correlogramas que no hay ningún rezago significativo (aparte del 0 del ACF, que es 1 por definición) que denote ningún tipo de estructura, por lo tanto podemos decir que los residuos son ruido blanco. Ahora, vamos a pintar también una serie de gráficas sobre los residuos:

ggtsdiag(arima6)

bp <- Box.test(arima6$residuals) # Test de Box-Pierce
bp$p.value

## [1] 0.6640325

lb <- Box.test(arima6$residuals, type="Ljung-Box") # Test de Ljung-Box
lb$p.value

## [1] 0.6610922

jb <- jarque.bera.test(arima6$residuals) # Test de Jarque-Bera
jb$p.value

##  X-squared 
## 0.04008699

sht<-shapiro.test(arima6$residuals) $ # Test de Shapiro-Wilk
sht$p.value

En R existe una función llamada auto.arima que puede calcular automáticamente cuál es la mejor combinación de órdenes para el modelo ARIMA:

auto.arima(co2ts, stepwise = FALSE, approximation = FALSE)

## Series: co2ts 
## ARIMA(1,0,1)(0,1,1)[12] with drift         
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.8999  -0.4391  -0.7686  0.1193
## s.e.  0.0479   0.0885   0.0930  0.0034
## 
## sigma^2 estimated as 0.09903:  log likelihood=-35.58
## AIC=81.16   AICc=81.58   BIC=96.18

En este caso la función sugiere que el mejor modelo que representa a la serie sería un ARIMA(1,0,1)(0,1,1) con tendencia, y que arroja un AIC de 81.16, siendo algo mejor al modelo que anteriormente habíamos visto, 83.83

Predicción de la serie

forecast1<-forecast(arima6, level = c(95), h = 50)
autoplot(forecast1)

Series Temporales