Introducción

Durante el documento se busca analizar y ajustar una serie de tiempo con componente estacional haciendo uso de la metodología Box y Jenkins (1976). La variable que contiene la información será llamada m3.

Metodología Box y Jenkins

Visualización de la serie de tiempo.

library(readr)
m3 <- read_csv("m3.cvs")
par(mfrow=c(2,1))
plot(m3$x,type = "l",main = "a) m3",ylab = "")
plot(m3$x[1:100],type = "o", main = "b) m3 (primeros 100 datos)",ylab = "")
Figura 1:Gráfico de la serie temporal

Figura 1:Gráfico de la serie temporal

La serie de tiempo está compuesta por 1000 datos, de los cuales en la Figura 1.a se puede observar una tendencia decreciente, mientras que la Figura 1.b que corresponden a los primeros 100 datos, pareciera presentar una componente estacional que se repite cada 7 periodos. Por otro lado, el comportamiento de dicha serie no muestra cambios en la varianza a lo largo del tiempo, con lo cual de golpe se podría suponer que su varianza es constante y que su media si se ve influenciada por el tiempo.

par(mfrow=c(2,1))
acf(m3$x,lag.max=35)
pacf(m3$x,lag.max = 35)
Figura 2:Gráfico ACF y PACF de m3

Figura 2:Gráfico ACF y PACF de m3

Para eliminar dichas tendencias como también se puede notar en el ACF de m3 con su decrecimiento lento (ver Figura 2), se plantean dos diferenciaciones:

  • La primera corresponde a la componente estacional con un lag \(1s\) (\(s=7\)) periodos.
  • La segunda a la parte no estacional con una lag de 1.

Con las acciones mencionadas anteriormente se generan los siguientes resultados:

x1_dif=diff(m3$x,lag = 1)
x2_dif=diff(m3$x,lag = 7)
x3_dif=diff(x2_dif,lag = 1)
par(mfrow=c(3,1))
plot(x1_dif,type = "l",ylab = "",main = "m3 con lag d=1")
plot(x2_dif,type = "l",ylab = "",main = "m3 con lag D=1")
plot(x3_dif,type = "l",ylab = "",main = "m3 con lag d=1 y D=1")
Figura 3:Gráfico de la serie temporal con los diferentes lags

Figura 3:Gráfico de la serie temporal con los diferentes lags

Los datos que han tenido las dos diferenciaciones (\(d=1\) y \(D=1\)) parecieran tener un comportamiento estacionario, por lo cual se partirá de estos para intentar ajustar un modelo.

ACF y PACF muestral

A continuación se ilustra la ACF y PACF muestral para proponer el modelo que mejor se ajuste a la serie de tiempo diferenciada, basado en lo que se logre identificar en ellas.

par(mfrow=c(2,1))
acf(x3_dif,col=rep(c(1,1,1,1,1,1,2),5),lag.max=35)
pacf(x3_dif,col=rep(c(1,1,1,1,1,1,2),5),lag.max = 35)
Figura 4:Gráfico ACF y PACF

Figura 4:Gráfico ACF y PACF

Se define como p, el mayor orden en el polinomio autorregresivo, y q es el mayor orden en el polinomio de media móviles.

Primero se analizan las líneas rojas de la Figura 4, las cuales darán indicios del comportamiento de la serie en la componente estacional, donde se puede observar que la PACF decrece y de forma sinusoidal, mientras que la ACF se anula después del primer rezago. Por tanto se pensaría que la componente estacional tiene un comportamiento ARIMA(0,1,1).

Dentro de cada periodo se identifica un comportamiento, en el ACF decrece lentamente y en el PACF sus lags son significativos solo en el primero, por lo cual se sospecha de un modelo ARIMA(1,1,0).

En definitiva se considera el modelo \(SARIMA(1,1,0)(0,1,1)_{7}\), del que se obtienen los siguientes resultados:

Diagnóstico del modelo

Para verificar la bondad del modelo, se utilizarán varios métodos con el fin de validar que se cumplan los supuestos. Como primera medida se calculan los residuales estandarizados:

sarima2<- arima(m3$x, order=c(1,1,0), seasonal=list(order=c(0,1,1),period=7),method = "CSS")
res_sarima2=residuals(sarima2)
res_est_sarima2=res_sarima2/(sarima2$sigma2^.5)
tsdiag(sarima2)
Figura 6: Gráficos de residuales

Figura 6: Gráficos de residuales 

qqnorm(res_est_sarima2,  xlab = "Cuantiles Teóreticos", ylab = "Cuantiles Muestrales")
qqline(res_est_sarima2,probs = c(0.25, 0.75))
Figura 7: QQnorm

Figura 7: QQnorm

Identificación de datos atípicos

Como se ilustra en la figura 6, no se evidencia un dato atípico en los residuales estandarizados del modelo, por lo cual procedemos a probar si es ruido blanco.

Normalidad

En la Figura 7 se puede percibir que se cumple el supuesto de normalidad al comparar los cuantiles teóricos con los cuantiles muestrales. Para robustecer dicha afirmación de emplea la prueba de Shapiro-Wilk.

shapiro.test(res_est_sarima2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_est_sarima2
## W = 0.99851, p-value = 0.5551

Que tiene el siguiente contraste de hipótesis, \(H_{0}:\) los residuales se distribuyen normal y \(H_{1}:\) los residuales no se distribuyen normal. Con el resultado de la prueba, no tenemos evidencia suficiente para rechazar el supuesto de normalidad y complementa lo visto en el QQplot.

Independencia

En la Figura 6, se pueden observar los valores de la prueba Ljung-Box la cual tiene el siguiente contraste \(H_{0}:\) los residuales se distribuyen de forma independiente \(H_{1}:\) los residuales no son independientes. De estos valores p se puede observar que todos están por arriba del nivel de significancia, lo que nos lleva a no rechazar la hipótesis nula, por tanto se cumple el supuesto, además se complementa esta afirmación viendo que el ACF de los residuales no tienen valores significativos.

Prueba de raíz unitaria

library(UnitCircle)              # contiene una función para obtener las raíces de los polinomios
x11()
uc.check(c(1, -coef(sarima2)[1]))
##        real complex outside
## 1 -1.217859       0    TRUE
## *Results are rounded to 6 digits.
Figura 8: Raíces unitarias

Figura 8: Raíces unitarias

Conclusiones

El modelo SARIMA propuesto se ajusta bien a los datos observados y se cumplen los supuestos. De modo que este modelo puede ser utilizado para conocer el comportamiento de nuestra serie temporal y para realizar pronósticos en el caso que se requiera.