INTRODUCCIÓN A LAS SERIES DE TIEMPO

Una serie tiempo es una secuencia de observaciones, medidos en determinados momentos del tiempo, ordenados cronológicamente y, espaciados entre sí de manera uniforme, así los datos usualmente son dependientes entre sí. El objetivo de una serie de tiempo es crear pronósticos a futuro de alguna variable \(X_t\) en la que \(t = 1, 2, 3, \cdots, n\) observaciones.

Componentes de una serie temporal

  • El análisis clásico de las series temporales se basa en la suposición de que los valores que toma la variable de observación es la consecuencia de tres componentes, cuya actuación conjunta da como resultado los valores medidos, estos componentes son:

  1. Componente tendencia: Se puede definir como un cambio a largo plazo que se produce en la relación al nivel medio, o el cambio a largo plazo de la media. La tendencia se identifica con un movimiento suave de la serie a largo plazo.

  2. Componente estacional: Muchas series temporales presentan cierta periodicidad o dicho de otro modo, variación de cierto período (semestral, mensual, etc.). Por ejemplo las Ventas al Detalle en Puerto Rico aumentan por los meses de noviembre y diciembre por las festividades navideñas. Estos efectos son fáciles de entender y se pueden medir explícitamente o incluso se pueden eliminar de la serie de datos, a este proceso se le llama desestacionalización de la serie.

  3. Componente aleatoria: Esta componente no responde a ningún patrón de comportamiento, sino que es el resultado de factores fortuitos o aleatorios que inciden de forma aislada en una serie de tiempo.

De estos tres componentes los dos primeros son componentes determinísticos, mientras que la última es aleatoria. Así se puede denotar la serie de tiempo como: \[X_t = T_t + E_t + I_t\] donde:

  • \(T_t\) es la tendencia.

  • \(E_t\) es la componente estacional.

  • \(I_t\) es la componente aleatoria.

Cargamos nuestras librerías a utilizar

library(lmtest)
library(readxl)
library(openxlsx)
library(forecast)
library(tseries)
library(timeSeries)
library(ggplot2)

Generamos nuestra base de datos

GDP = read_xls("C:\\Users\\Brook\\Desktop\\GDP.xls")

Creamos nuestra variable

Como nuestros datos se tratan del principal índice económico, será fácil de utilizar y lograr observar con detenimiento. Generaremos nuestros datos desde la primera observación obtenida hasta la última, esto también lo haremos de manera trimestral como se observa a continuación:

GDP = ts(GDP, start = c(1947, 1), frequency = 4)
autoplot(GDP, color = "darkgoldenrod4", size = 0.7) +
  labs(title = "Producto Interno Bruto (1947/T1 - 2020/T2)",
       subtitle = "Estados Unidos, en Millones de USD (Precios Corrientes)", caption = "FRED (2021) https://fred.stlouisfed.org/series/GDP") +
  xlab("Años") + ylab("Producto Interno Bruto")

Una de las primeras cuestiones de las series de tiempo es que nuestras variables deben de ser estacionarias, en otras palabras, que la media y la varianza sean constantes y no afecte el tiempo dentro de ellas:

  • Propiedades a cumplir:

\[Media: E(X_t)=E(X_t+_k)=\mu\] \[Varianza:V(X_t)=V(X_t+_k)=\sigma^2\] \[Covarianza:\gamma_k=E[(X_t-\mu)(X(_t+_k)-\mu)]\]

Serie estacionaria y con estacionalidad

Ahora lo que haremos será crear el primer diferencial de la serie original para que así puedamos trabajar con el pronóstico a futuro y sin ningun problema de raíz unitaria. También, trabajaremos con logaritmo para que así pasemos de niveles a tasas la serie:

logGDP1 = diff(log(GDP))

autoplot(logGDP1, color = "darkgoldenrod4", size = 0.7) + 
  labs(title = "Primer Diferencial del Logaritmo del PIB (1947/T1 - 2020/T4)", subtitle = "Estados Unidos (Tasas)", caption = "FRED (2021)") +
  xlab("Años") + ylab("Diferencial")

Observamos que la serie ahora sí se ha vuelto estacionaria y estacional al generar su primer diferencial del logaritmo. Debemos de comprobarlo a través de las pruebas de Dickey-Fuller aumentada y Phillips-Perron para comprobar que realmente lo sea.

Revisión de estacionalidad

adf.test(logGDP1)
## Warning in adf.test(logGDP1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  logGDP1
## Dickey-Fuller = -5.3928, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
pp.test(logGDP1)
## Warning in pp.test(logGDP1): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  logGDP1
## Dickey-Fuller Z(alpha) = -252.24, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

En ambos casos observamos que el valor p de Pearson es menor que el valor crítico de 0.05, por lo que podemos concluir que rechazamos la Hipótesis nula. En otras palabras, la serie es estacionaria y con estacionalidad y no existe una raíz unitaria que afecte a nuestro modelo. El tiempo no influye dentro de nuestra serie de tiempo.

Autocorrelación y Autocorrelación parcial

Autocorrelación En ocasiones en una serie de tiempo acontece, que los valores que toma una variable en el tiempo no son independientes entre sí, sino que un valor determinado depende de los valores anteriores, existen dos formas de medir esta dependencia de las variables

Autocorrelación Parcial La autocorrelación parcial mide la correlación entre dos variables separadas por \(k\) periodos cuando no se considera la dependencia creada por los retardos intermedios existentes entre ambas.

autoplot(acf(logGDP1), color = "deepskyblue", size = 1.5) + 
  labs(title = "Autocorrelación",
       subtitle = "95% de confianza") +
  xlab("Rezagos") +
  ylab("Autocorrelación")

autoplot(pacf(logGDP1), color = "deepskyblue", size = 1.5) +
  labs(title = "Autocorrelación Parcial",
       subtitle = "95% de confianza") +
  xlab("Rezagos") +
  ylab("Autocorrelación Parcial")

Descomposición de la serie de tiempo

Ahora descompondremos la serie para obtener sus componentes. Como lo es: * La serie origninal * Tendencia * Estacionalidad * Recordatorio

autoplot(decompose(logGDP1)) + labs(title = "Descomposición de las adiciones
de la serie de tiempo", 
subtitle = "(Datos) (Tendencia) (Estacionalidad) (Recordatorio)")

Generar el modelo ARIMA

El primer paso será generar nuestro modelo ARMA sin estacionalidad y también sin estacionariedad, esto para ver el mejor comportamiento que se tiene para la serie de tiempo y ver cuántas Autorregresiones (AR), Integraciones (I) y Medias Móviles (MA) se necesitan para el mejor pronóstico:

modelo1 = auto.arima(logGDP1, seasonal = FALSE, stationary = FALSE)

summary(modelo1)
## Series: logGDP1 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.1355  -0.9536
## s.e.  0.0659   0.0176
## 
## sigma^2 estimated as 0.000162:  log likelihood=862.84
## AIC=-1719.68   AICc=-1719.6   BIC=-1708.64
## 
## Training set error measures:
##                         ME       RMSE         MAE     MPE     MAPE      MASE
## Training set -0.0005473788 0.01266207 0.007428117 34.1741 111.2191 0.6840467
##                    ACF1
## Training set -0.0169306
autoplot(modelo1) + labs(title = "Círculos de Raíces Unitarias",
                         subtitle = "Raíces inversas del modelo AR MA",
                         caption = "Modelo sin estacionalidad y sin ser estacionaria") +
  xlab("Reales") + 
  ylab("Imaginarios")

checkresiduals(modelo1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 11.91, df = 6, p-value = 0.06401
## 
## Model df: 2.   Total lags used: 8

Observamos que el mejor modelo ARIMA es de orden (1,1,1), con los coeficientes \(AR=0.1355\) y \(MA=-0.9536\) a una integración más. Hay que recordar que ya habíamos integrado la serie una sola vez anteriormente para hacerla estacionaria, por lo que el ARMIA propuesto es de orden (1,2,1). De igual forma, observamos que en la prueba Ljung-Box muestra un p-value mayor al crítico de 0.06, por lo que aceptamos la Hipótesis nula, los residuales se distribuyen de manera normal y son estacionarios.

Pronóstico del ARIMA

Generamos entonces un pronóstico a 20 observaciones futuras el comportamiento de nuestra serie de tiempo de nuestro modelo ARIMA. Lo que viene siendo 5 años más ó 20 trimestres más.

pronostico1 = forecast(modelo1, h = 20)

autoplot(pronostico1) +
  labs(title = "Pronóstico: 20 trimestres para el Logaritmo del Primer Diferencial del PIB",
       subtitle = "ARIMA (1,2,1), Sin estacionalidad y sin estacionariedad") +
  xlab("Años") + ylab("Primer diferencial ")

Observamos que el comportamiento es algo lineal y no muestra grandes cambios a futuro con nuestra serie de tiempo del PIB de Estados Unidos. Esto puede deberse a que sea un modelo sin estacionalidad ni que sea estacionaria.

Generar el modelo SARIMA

El modelo SARIMA es una variación a lo que vienen siendo un modelo ARIMA, ya que estos modelos aparte de incluir la parte (p,d,q) incluye la parte (P,D,Q), lo estacional. A lo que el pronóstico a futuro es mejor como se observ a continuación:

modelo2 = auto.arima(logGDP1)

summary(modelo2)
## Series: logGDP1 
## ARIMA(1,1,3)(2,0,0)[4] 
## 
## Coefficients:
##          ar1      ma1     ma2      ma3     sar1     sar2
##       0.1107  -1.0114  0.2911  -0.2197  -0.1103  -0.1072
## s.e.  0.2833   0.2844  0.2765   0.0826   0.0887   0.0812
## 
## sigma^2 estimated as 0.0001584:  log likelihood=868.07
## AIC=-1722.14   AICc=-1721.75   BIC=-1696.38
## 
## Training set error measures:
##                         ME      RMSE         MAE      MPE     MAPE      MASE
## Training set -0.0005169623 0.0124331 0.007285186 30.78685 103.9186 0.6708844
##                      ACF1
## Training set -0.002714013
autoplot(modelo2) + labs(title = "Círculos de Raíces Unitarias",
                         subtitle = "Raíces inversas del modelo AR MA",
                         caption = "Modelo estacional y estacionario") +
  xlab("Reales") + 
  ylab("Imaginarios")

checkresiduals(modelo2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,3)(2,0,0)[4]
## Q* = 5.5734, df = 3, p-value = 0.1343
## 
## Model df: 6.   Total lags used: 9

Vemos que se trata de un modelo SARIMA (1,1,3)(2,0,0), por lo que nuestra ecuación de pronóstico está formada de la siguiente forma: \(AR_1=0.1107\), \(MA_1=-1.0114\), \(MA_2=0.2911\), \(MA_3=-0.2197\), \(sAR_1=-0.1103\), y \(sAR_2=-0.1072\). De igual forma observamos que la prueba Ljung-Box el p-value es mayor al crítico de 0.05, de esta forma aceptamos la Hipótesis nula y concluimos que los errores se distribuyen de manera normal y llegan a ser estacionarios.

Pronóstico

De igual forma, crearemos un pronóstico para el futuro a 20 trimestre en adelante:

pronostico2 = forecast(modelo2, h = 20)

autoplot(pronostico2) +
  labs(title = "Pronóstico: 20 trimestres para el Logaritmo del Primer Diferencial del PIB",
       subtitle = "SARIMA (1,1,2), (2,0,0), Con estacionalidad y estacionariedad") +
  xlab("Años") + ylab("Primer diferencial")

Vemos que el comportamiento predictorio es aun mejor en el SARIMA que en el ARIMA.

Bibliografía