Taller SARIMA

Librerías usadas

library(forecast)
library(ggplot2)
library(tseries)
library(Metrics)
library(car)
library(nortest)
library(lmtest)

Exploración de la serie

Analizaremos la serie de tiempo que contiene los registros mensuales de la producción de cemento en Colombia, a partir de enero de 2010 hasta diciembre de 2023.

Para empezar vamos a leer y a definir la serie.

datos <- read.table("C:/Users/danii/OneDrive/Escritorio/Simposio/series.txt",header=T)
datos <- datos$cemento1
serie <- ts(datos, frequency = 12, start = c(2010,1), end = c(2023,12))

Ahora, pensando en realizar validación cruzada, vamos a definir una parte de la serie como los datos de entrenamiento y los demás como datos de prueba. En este caso, la base de entrenamiento será la que va desde enero del 2010 hasta diciembre del 2022, por su parte, el último año será la base de prueba.

N = length(serie)
train = ts(serie, frequency=12, start=c(2010,1), end = c(2022,12)) 
test = ts(serie[(N-12+1):N], frequency=12, start=c(2023,1), end=c(2023,12)) 

Identificación de componentes

Lo primero que haremos será analizar la serie, de esta manera podremos ir pensando en qué tipo de modelos sería adecuado probar.

autoplot(train)

Lo primero que notamos es que la serie parece tener una tendencia al alza. La serie al inicio de los registros toma valores cercanos al 800000, luego se va evidenciando un paulatino incremento conforme pasan los años, para los últimos años la serie toma valores alrededor de 1200000. La media de la serie no es estacionaria, es decir, esta no permanece constante y parece en cambio depender del tiempo.

ggseasonplot(train, polar = T)

La serie parece tener un componente estacional, esto no es tan claro pero algunos detalles en el comportamiento de la serie, tales como la recurrencia de valores más bajos en el mes de enero y valores más uniformes en la temporada del mes de julio al mes de octubre, llevan a pensar que esto obedece a un comportamiento estacional.

plot(decompose(train, type = 'additive')) 

Por último, observamos que no es muy claro si la variación de la serie cambia con el tiempo. Aparte del gran pico descendente en el primer semestre del año 2020 (explicado por la pandemia del COVID-19), no es muy claro si los picos aumentan o disminuyen en su magnitud. El gráfico de la serie descompuesta nos muestra un patrón estacional más claro, una tendencia al crecimiento y el componente aleatorio que resta, este parece no variar, sin embargo, no es certero esta observación, por lo que en el futuro se tienen en cuenta las dos posibilidades (que la varianza sea o no estacionaria)

Recurrimos al análisis de la gráfica de la función de autocorrelación ACF y de la función de autocorrelaciones parciales PACF para corroborar nuestras hipótesis.

ggAcf(train, lag.max = 50)
ggPacf(train, lag.max = 50)

La gráfica de la ACF nos muestra que el comportamiento de sus rezagos es de decaimiento leve, esto es, los valores no caen rápidamente en los primeros rezagos dentro de las bandas y en cambio parece que sigue un comportamiento leve e irregular en forma de “U”. También es evidente que los rezagos estacionales (múltiplos de 12) son más altos que los demás rezagos.

Por otro lado, la PACF muestra un comportamiento aparentemente sinusoidal, pues los valores oscilan de arriba a abajo del cero. Los dos primeros regazos regulares se salen de las bandas y el primer rezazgo estacional también

Ajustes de modelos

Ahora, después de haber identificado el comportamiento que presenta la serie, procedemos con el ajuste de los modelos. En este caso se probarán tres tipos de modelos: los modelos Holt-Winters, los modelos Box-Jenkins y la red neuronal autorregresiva estacional. A continuación se presentan los modelos y la comparación entre sí.

Holt-Winters

Box-Jenkins

Transformación y diferenciación

Para poder pensar en un modelo Box-Jenkins debemos garantizar que la serie sea estacionaria. Esto se logra eliminando la tendencia y la estacionalidad. También se procura estabilizar la varianza de la serie a lo largo del tiempo, aplicando una transformación a la serie en los casos que sean necesarios y justificados.

Primero se va a evaluar si es necesario estabilizar la varianza de la serie por medio de una transformación, en este caso vamos a aplicar la transformación Box-Cox.

lambda <- BoxCox.lambda(train)
lambda

[1] -0.2737781

En este caso identificamos un valor de -0.2737781 para lambda, el parámetro de la transformación Box-Cox.

train_T <- BoxCox(train, lambda)
plot(train)
plot(train_T)

No hay un cambio evidiente entre la serie original y la serie transformada. Por tal razón, siguiendo la recomendación de algunos autores de no aplicar una transformación a menos que sea estrictamente necesaria, se decide trabajar con la serie original.

Continuamos con la diferenciación de la serie. Para series de tiempo que necesiten de diferenciación estacional y regular (series con tendencia y estacionalidad), se pueden aplicar en el orden que el analista lo decida. En este caso se decide aplicar una diferenciación regular primero.

Para empezar, identificamos el número de diferenciaciones estacionales que debemos aplicar.

ndiffs(train)

[1] 1

Procedemos a diferenciar la serie con respecto a la longitud del rezago estacional.

# data_T_DS <- diff(data_T, lag = frequency(data))
train_d = diff(train)
autoplot(train_d)

La serie a simple vista parece ser ya estacionaria, sin embargo, no sabemos si la diferenciación regular eliminó el componente estacional que habíamos identificado. Procedemos entonces a analizar las gráficas de ACF y PACF de la serie diferenciada.

ggAcf(train_d, lag.max = 50)
ggPacf(train_d, lag.max = 50)

La prueba ACF deja ver que la tendencia fue eliminada, entonces no hay necesidad de otra diferenciación regular. Por otro lado, también muestra claramente unos picos más altos en los rezagos estacionales, lo que nos confirma que la estacionalidad no ha sido eliminada y que es necesario aplicar una diferenciación estacional. A continuación aplicamos la diferenciación estacional.

train_d_s <- diff(train_d, lag = 12)

Realizada las diferenciaciones, para confirmar que la serie no es estacionaria y que se puede proceder a pensar en el orden del modelo SARIMA aplicamos la prueba de la raíz unitaria Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

Recordemos que la prueba KPSS presenta las siguientes hipótesis:

\[ H_0 = La\;serie\;es\;estacionaria \] \[ H_1 = La\;serie\;no\;es\;estacionaria \]

kpss.test(train_d_s)

Warning in kpss.test(train_d_s): p-value greater than printed p-value

    KPSS Test for Level Stationarity

data:  train_d_s
KPSS Level = 0.022813, Truncation lag parameter = 4, p-value = 0.1

El valor p de la prueba KPSS tiene un valor más grande que 0.1, por lo tanto, aceptamos la hipótesis nula que dice que la serie es estacionaria. Por medio de estas dos diferenciaciones se observa un proceso estacionario a partir del cuál podemos empezar a concluir sobre los modelos SARIMA.

Identificación de parámetros

Como estamos trabajando con una serie con componente estacional, debemos pensar en los parámetros p, q, P y Q. Para esto graficamos la ACF y la PACF de la serie estacionaria.

ggAcf(train_d_s, lag.max = 50)
ggPacf(train_d_s, lag.max = 50)

Para determinar los parámetros P y Q analizamos los rezagos estacionales, estos son los múltiplos de la frecuencia de la serie, es decir 12, 24, 36 y 48.
Para el caso de la ACF, se ve que solo el primer rezago estacional se sale de las bandas, en cambio en el PACF los rezagos estacionales tienen un decaimiento lento. Esto nos lleva a considerar un posible MA(1) estacional.

Continuando con los patrones regulares, estudiamos los rezagos que hacen parte de un periodo estacional. De estos se puede notar que en la ACF solo el primer rezago sale de las bandas y la PACF tiene un comportamiento sinusoidal. En consecuencia, podemos pensar en el modelo MA(1).

De este análisis y teniendo en cuenta las diferenciaciones hechas, el modelo que se propone es un SARIMA(0,1,1)(0,1,1)

Identificación modelos

Ahora, habiendo previamente identificado los componentes de manera manual. Se usa la función auto.arima para identificar otro posible modelo. Esto lo haremos sobre la serie sin diferenciar y usaremos el criterio BIC.

auto.arima(train,  ic=c("bic"), approximation=FALSE)

Series: train 
ARIMA(0,1,1)(1,0,0)[12] 

Coefficients:
          ma1    sar1
      -0.6497  0.4765
s.e.   0.0758  0.0703

sigma^2 = 5.209e+09:  log likelihood = -1954.7
AIC=3915.4   AICc=3915.56   BIC=3924.53

La función auto.arima identificó el modelo SARIMA(0,1,1)(1,0,0). Se comparan los dos modelos SARIMA identificados con el objetivo de dejar el mejor y compararlo con los otros métodos propuestos.

• Modelo 1: SARIMA(0,1,1)(0,1,1)

• Modelo 2: SARIMA(0,1,1)(1,0,0)

Ajuste de los modelos

Vamos a estimar los coeficientes para cada modelo y a probar su significancia usando la prueba de la normal.

Recordemos que esta prueba presenta las siguientes hipótesis:

Para el caso de los coeficientes \(\phi_i\)

\[ H_0: \phi_i=0 \]

\[ H_1: \phi_i\neq 0 \]

Para el caso de los coeficientes \(\theta_j\)

\[ H_0: \theta_j=0 \]

\[ H_1: \theta_j\neq 0 \]

Para los coeficientes estacionales es igual, la diferencia es que en el subíndice del órden también se agrega una “s”

• Modelo 1: SARIMA(0,1,1)(0,1,1)

sarima1 <- Arima(train, order=c(0,1,1), seasonal=c(0,1,1))
summary(sarima1)

Series: train 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.5838  -0.8395
s.e.   0.0883   0.0898

sigma^2 = 4.494e+09:  log likelihood = -1798.53
AIC=3603.06   AICc=3603.23   BIC=3611.95

Training set error measures:
                    ME     RMSE      MAE        MPE     MAPE    MASE       ACF1
Training set -2093.309 63735.36 42280.56 -0.5149275 4.182649 0.56986 0.03966277

coeftest(sarima1)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.583815   0.088295 -6.6121 3.789e-11 ***
sma1 -0.839513   0.089792 -9.3495 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analizamos el valor p que tiene cada coeficiente ajustado y concluimos la validez de aquellos coeficientes donde el valor p sea menor que la significancia.

• \(\theta_1\) (ma1) ⟶ Como el valor p es menor que la significancia \(\alpha = 0.05\), se concluye que es significativo.

• \(\theta_{1s}\) (sma1) ⟶ Como el valor p es menor que la significancia \(\alpha = 0.05\), se concluye que es significativo.

De este modelo se concluye que todos los parámetros son significativos. Por lo tanto, el modelo también lo es.

• Modelo 2: SARIMA(0,1,1)(1,0,0)

sarima2 <- Arima(train, order=c(0,1,1), seasonal=c(1,0,0))
summary(sarima2)

Series: train 
ARIMA(0,1,1)(1,0,0)[12] 

Coefficients:
          ma1    sar1
      -0.6497  0.4765
s.e.   0.0758  0.0703

sigma^2 = 5.209e+09:  log likelihood = -1954.7
AIC=3915.4   AICc=3915.56   BIC=3924.53

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
Training set 5420.234 71474.34 52255.73 0.1510234 5.228262 0.7043059 0.04881163

coeftest(sarima2)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.649714   0.075759 -8.5761 < 2.2e-16 ***
sar1  0.476514   0.070279  6.7803 1.199e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analizamos el valor p que tiene cada coeficiente ajustado y concluimos la validez de aquellos coeficientes donde el valor p sea menor que la significancia.

• \(\theta_1\) (ma1) ⟶ Como el valor p es menor que la significancia \(\alpha = 0.05\), se concluye que es significativo.

• \(\phi_{1s}\) (sar1) ⟶ Como el valor p es menor que la significancia \(\alpha = 0.05\), se concluye que es significativo.

De este modelo se concluye que todos los parámetros son significativos. Por lo tanto, el modelo también lo es.

Validación de supuestos

Para las pruebas estadísticas se tendrá en cuenta un nivel de significancia \(\alpha=0.01\)

checkresiduals(sarima1, test = FALSE)

checkresiduals(sarima2, test = FALSE)

Los residuales de ambos modelos parecen ruido blanco. Aunque si comparamos las gráficas de la ACF, miramos que los residuales del segundo modelo parecen estar mejor que los del primero, ya que ningún rezago está fuera de las bandas, contrario a los residuales del primer modelo donde algunos rezagos (especialmente los estacionales) se salen.

Supuesto de normalidad

Modelo 1

gghistogram(sarima1$residuals)
boxplot(sarima1$residuals)
qqPlot(sarima1$residuals)

[1] 123 137

Las pruebas gráficas nos muestran que, a pesar de que en principio parece haber una normalidad, en el qqplot se sugiere un problema de colas pesadas, pues estas se salen de los valores teóricos de normalidad.

\[ H_0:Los\;residuales\;siguen\;una\;distribución\;normal \]

\[ H_1:Los\;residuales\;no\;siguen\;una\;distribución\;normal \]

shapiro.test(sarima1$residuals)

    Shapiro-Wilk normality test

data:  sarima1$residuals
W = 0.89512, p-value = 4.218e-09

ad.test(sarima1$residuals)

    Anderson-Darling normality test

data:  sarima1$residuals
A = 3.4336, p-value = 1.268e-08

cvm.test(sarima1$residuals)

    Cramer-von Mises normality test

data:  sarima1$residuals
W = 0.54292, p-value = 1.086e-06

Los residuales del modelo no pasan ninguna prueba de normalidad.

Modelo 2

gghistogram(sarima2$residuals) 
boxplot(sarima2$residuals) 
qqPlot(sarima2$residuals)

[1] 123 125

Las pruebas gráficas nos muestran que, a pesar de que en principio parece haber una normalidad, en el qqplot se sugiere un problema de colas pesadas, pues estas se salen de los valores teóricos de normalidad. Algo notable es que el boxplot refleja menos outliers que aquellos del modelo 1.

\[ H_0:Los\;residuales\;siguen\;una\;distribución\;normal \]

\[ H_1:Los\;residuales\;no\;siguen\;una\;distribución\;normal \]

shapiro.test(sarima2$residuals)

    Shapiro-Wilk normality test

data:  sarima2$residuals
W = 0.9695, p-value = 0.001569

ad.test(sarima2$residuals)

    Anderson-Darling normality test

data:  sarima2$residuals
A = 1.0623, p-value = 0.008402

cvm.test(sarima2$residuals)

    Cramer-von Mises normality test

data:  sarima2$residuals
W = 0.17802, p-value = 0.0101

A pesar de que el valor de las pruebas aplicadas mejora bastante con respecto a aquellos del primer modelo, estos residuales tampoco pasan ninguna prueba de normalidad.

Supuesto de homocedasticidad

Modelo 1

df <- data.frame(Fitted = as.numeric(fitted(sarima1)),
                 Residuals = as.numeric(residuals(sarima1)))

ggplot(df, aes(x = seq_along(Residuals), y = Residuals)) +
  geom_point() +
  labs(x = "Índice", y = "Residual")

ggplot(df, aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Valores Ajustados", y = "Residuales") 

Gráficamente los residuales parecen tener varianza no constante.

\[ H_0:La\;varianza\;de\;los\;residuales\;es\;constante\; \]

\[ H_1:La\;varianza\;de\;los\;residuales\;no\;es\;constante\; \]

bptest(residuals(sarima1)~fitted(sarima1))

    studentized Breusch-Pagan test

data:  residuals(sarima1) ~ fitted(sarima1)
BP = 3.2504, df = 1, p-value = 0.07141

El modelo uno pasa la pueba Breusch-Pagan, por lo que se cumple el supuesto de homocedasticidad.

Modelo 2

df <- data.frame(Fitted = as.numeric(fitted(sarima2)),
                 Residuals = as.numeric(residuals(sarima2)))


ggplot(df, aes(x = seq_along(Residuals), y = Residuals)) +
  geom_point() +
  labs(x = "Índice", y = "Residual")
ggplot(df, aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Valores Ajustados", y = "Residuales") 

Gráficamente los residuales parecen tener varianza constante.

\[ H_0:La\;varianza\;de\;los\;residuales\;es\;constante\; \]

\[ H_1:La\;varianza\;de\;los\;residuales\;no\;es\;constante\; \]

bptest(sarima2$residual~sarima2$fitted)

    studentized Breusch-Pagan test

data:  sarima2$residual ~ sarima2$fitted
BP = 0.55121, df = 1, p-value = 0.4578

Se concluye que la varianza de los residuales es constante para el modelo 2.

Supuesto de independencia

Modelo 1

ggAcf(sarima1$residuals)
ggPacf(sarima1$residuals)

La gráfica muestra una aparente independencia entre los residuales, aunque algunos rezagos de la ACF y PACF salen de las bandas.

\[ H_0:Los\;residuales\;son\;independientes\; \]

\[ H_1:Los\;residuales\;no\;son\;independientes \]

Box.test(x = residuals(sarima1), lag=36, type="Ljung-Box")

    Box-Ljung test

data:  residuals(sarima1)
X-squared = 36.771, df = 36, p-value = 0.433

El valor p es mayor que el nivel de significancia. Se acepta la hipótesis nula y se valida el supuesto de independencia de los residuales.

Modelo 2

ggAcf(sarima2$residuals)
ggPacf(sarima2$residuals)

La gráfica muestra una aparente independencia entre los residuales, no hay ningún rezago que se salga de las bandas.

\[ H_0:Los\;residuales\;son\;independientes\; \]

\[ H_1:Los\;residuales\;no\;son\;independientes \]

Box.test(x = residuals(sarima2), lag=36, type="Ljung-Box")

    Box-Ljung test

data:  residuals(sarima2)
X-squared = 30.071, df = 36, p-value = 0.7458

El valor p es mayor que el nivel de significancia. Se acepta la hipótesis nula y se valida el supuesto de independencia de los residuales.

Comparación modelos

Usaremos como medidas de desempeño el Mean Absolute Percentual Error (MAPE), el criterio de información de Akaike AIC y el criterio bayesiano de Schwarz BIC.

MAPEs = cbind(MAPEmod1 = mape(sarima1$x, sarima1$fitted),
              MAPEmod2 = mape(sarima2$x, sarima2$fitted))
MAPEs

       MAPEmod1   MAPEmod2
[1,] 0.04182649 0.05228262

AICs = cbind(AICmod1 = AIC(sarima1), 
             AICmod2 = AIC(sarima2)) 
AICs

      AICmod1 AICmod2
[1,] 3603.059  3915.4

BICs = cbind(BICmod1 = BIC(sarima1),
             BICmod2 = BIC(sarima2)) 
BICs

      BICmod1  BICmod2
[1,] 3611.947 3924.531

Por unanimidad el mejor modelo es el modelo 1, pues para todas las medidas de desempeño tenidas en cuenta es el que obtuvo un valor más bajo. Esto nos asegura que sus predicciones son más confiables.

Red neuronal