Proyecto Series de Tiempo

Analisis Descriptivo

Realizamos un gráfico para la serie de tiempo el cual nos muestra el rendimiento anual de grano en el campo Broadbalk en Rothamsted a lo largo del periodo 1852-1925. Con la sentencia sum(is.na(ts_38)) detectamos que no hay NA’s en nuestra serie.

Serie de tiempo

Se aprecian bastantes datos atípicos, algunos en los años 1863, 1876, 1883, 1904 y 1920, por mencionar algunos. Estos datos son normales por las diversas eventualidades presentes en el sector agrícola que provocan una baja de gran impacto en el rendimiento del grano (fenómenos meteorológicos, plagas, descuido del cultivo, etc.) o, por otro lado, un alza (climas estacionales ideales, implementación de métodos de cultivo eficientes, etc.)

Datos extremos en la serie de tiempo

La tendencia del rendimiento del grano tiene un comportamiento variable, siendo alcista en algunos periodos y bajista en algunos otros.
A simple vista, resulta complicado ver si la serie tiene algún efecto estacional, pero al realizar los gráficos de los correlogramas, se logra apreciar que la serie no presenta dicho efecto pues se observa un comportamiento distinto durante todo el periodo de tiempo.

Gráficos de autocorrelación

Finalmente, la varianza indica ser constante a lo largo del tiempo, pero no podemos omitir las enormes variaciones originadas por la gran cantidad de datos atípicos, así como la significativa baja que hay en los datos entre 1860 y 1875.
Por todas las observaciones anteriores, concluimos que la serie es no estacionaria por lo cual habrá que realizar ciertas transformaciones y un proceso de diferenciación para que pueda serlo según los requerimientos del desarrollo del presente trabajo.

Descomposición Clásica

Para comenzar con la descomposición clásica de la serie de tiempo, recordemos que no contamos con un componente estacional, por lo que podría resultar más significativo considerar que la serie sigue un modelo aditivo. Sin embargo, realizar una afirmación al respecto resulta complicado, pues no podemos omitir lo observado en el análisis previo sobre la varianza. Así que decidimos usar la transformación LN para reducir la varianza de manera general y así asegurarnos de que la serie sigue un modelo aditivo. 

Tendencia

Aplicamos el método clásico “Suavizamiento a través de Medias Móviles Simples” para obtener la componente de tendencia. Utilizamos la función nativa de R “ma(x, order)” la cual sustituye cada valor de la serie por la media obtenida de los datos dentro de un intervalo centrado en la observación a sustituir según el valor order que coloquemos (influyendo si es par o impar), en este caso consideramos que order = 8 es un valor adecuado para obtener promedios con suficiente información. Al final de este proceso, es importante hacer constantes los valores que no fueron sustituidos, utilizando el primer y último dato, según sea el caso.

Tendencia

Con esta información es más claro ver las fluctuaciones que tiene la tendencia, confirmando que de positiva pasa a negativa y viceversa aunque con distintas magnitudes según el periodo que observemos. De igual forma, notamos que incluso después de suavizar, siguen persistiendo algunas pequeñas variaciones, sin embargo, no son de igual relevancia como en el gráfico original de la serie.

Change Points

A pesar de que la tendencia cambia a lo largo del tiempo, podemos asegurar que no hay Change Points que deriven en algún cambio abrupto en ella y provoquen la necesidad de seccionar nuestra serie en partes para poder realizar predicciones de manera coherente al tomar toda la serie como base.

Change Points

Componente aleatorio

Ahora sólo nos falta obtener el componente aleatorio de la serie, el cual obtenemos restando el componente de tendencia a la serie transformada. Finalmente tomamos, de manera ilustrativa, presentamos una muestra de los valores obtenidos de la descomposición, así como las respectivas gráficas.

## Time Series:
## Start = 1852 
## End = 1856 
## Frequency = 1 
##       LOGts_38 tendencia    aleatorio
## 1852 0.0000000 0.5324359 -0.532435888
## 1853 0.2851789 0.5324359 -0.247256946
## 1854 0.5306283 0.5324359 -0.001807637
## 1855 0.8796267 0.5324359  0.347190859
## 1856 0.2390169 0.5324359 -0.293418988

Descomposición clásica de la serie de tiempo.

Ajuste de Series

Para realizar el ajuste del modelo, aplicaremos una metodología de validacion cruzada por lo cual, en primer lugar, separamos la serie de tiempo en dos conjuntos, el de entrenamiento (training) y el de evaluación (testing) para analizar el primero y asi emplear un modelo. Utilizamos una proporción 80%, 20% para los conjuntos anteriores.

Dentro de la metodología aprendida para encontrar un modelo adecuado, graficamos los correlogramas.

Correlogramas

Observamos que la estructura de correlación se preserva mayormente a la de la serie completa por lo que consideramos que debido a esa fuerte correlación entre lag’s adyacentes, la varianza no es constante a través del tiempo, la tendencia es deterministica y existen outliers significativos; inferimos que la serie es no estacionaria.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_training
## Dickey-Fuller = -2.6744, Lag order = 3, p-value = 0.3021
## alternative hypothesis: stationary

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

## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_training
## KPSS Level = 0.14017, Truncation lag parameter = 3, p-value = 0.1

Ante la situacion, buscamos reducir el impacto de los outliers, reducir la varianza y hacer lo más constante posible la tendencia de modo que la serie resulte lo más cercano a una serie debilmente estacionaria.

Transformación

Para este trabajo empleamos 3 distintas transformaciones con el objetivo de lograr el ajuste más eficiente, las cuales son: BOX-COX, LOG y SQRT. Cada una logra distintos objetivos, la primera es la más neutra en cuanto a varianza constante, tendencia constante y reducir los shocks de la serie (reduce outliers), la segunda vuelve constante la varianza salvo un intervalo puntual y logra hacer constante la tendencia (casi cero), la tercera es la que más reduce la varianza total y el outlier no es tan grande como en la transformación logaritmo.

Transformaciones

A continuación verificamos la estacionareidad de las series transformadas bajo dos distintos enfoques.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(ts_training)
## Dickey-Fuller = -2.8425, Lag order = 3, p-value = 0.2341
## alternative hypothesis: stationary

## Warning in kpss.test(log(ts_training)): p-value greater than printed p-value

## 
##  KPSS Test for Level Stationarity
## 
## data:  log(ts_training)
## KPSS Level = 0.21521, Truncation lag parameter = 3, p-value = 0.1

De lo anterior, se puede ver que la serie transformada con LOG no es estacionaria pues no rechazamos en Augmented Dickey-Fuller y rechazamos en la KPSS, es decir, hay evidencia a favor de la existencia de al menos una raíz unitaria y no es trend-stationary. De la misma forma concluimos con las transformaciones BOX-COX y SQRT.

Para conseguir una serie estacionaria, buscamos la(s) diferencia(s) adecuada(s) mediante el metodo de varianzas expuesto en clase.

Diferenciamos para quitar la tendencia y elegimos la de menor varianza.

  var(diff(log(ts_training)))

## [1] 0.5879354

  var(diff(log(ts_training), 2))

## [1] 0.2510723

  var(diff(diff(log(ts_training))))

## [1] 2.060346

  var(diff(diff(log(ts_training)),2))

## [1] 0.4783001

  #Campeona
  var(diff(log(ts_training), 2))

## [1] 0.2510723

Ahora, aplicando la diferenciación ganadora, probamos estacionariedad.

## Warning in adf.test(diff(log(ts_training), 2)): p-value smaller than printed p-
## value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(log(ts_training), 2)
## Dickey-Fuller = -5.5141, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

## Warning in kpss.test(diff(log(ts_training), 2)): p-value greater than printed p-
## value

## 
##  KPSS Test for Level Stationarity
## 
## data:  diff(log(ts_training), 2)
## KPSS Level = 0.090398, Truncation lag parameter = 3, p-value = 0.1

Como pasamos la primera prueba de estacionareidad, Augmented Dickey-Fuller, que es la más potente, aunque se rechace en la KPSS, consideramos que la anterior es una serie estacionaria.

Ahora ya con la transformación adecuada, nos encaminamos en la busqueda del modelo, la diferencia aplicada nos indica que en el sARIMA de ajuste, convendría un parámetro de D=1 con s=2.

Para visualizar la correlación de esta nueva serie y así hallar los parametros faltantes, graficamos los correlogramas.

Autocorrelogramas de la transformación diferenciada

Los correlogramas resultan difusos, pero observamos que en el ACF Y PACF, en el lag 2 se cortan bruscamente, y en el 4 también, no hay un claro decaimiento exponencial o senoidal, por lo que consideraremos agregar ambos terminos SAR,AR,MA,SMA en un principio.

Tras comparacion exhaustiva de aproximadamente 25 modelos (considerando transformaciones, diferencias) en la cual tomamos en cuenta el cumplimiento de los supuestos y principalmente RMSE Y MAE del ajuste en el conjunto training, RMSE Y MAE del ajuste/forecast en el conjunto testing; y los criterios de informacion AIC, AICc y BIC), obtuvimos el mejor modelo en general, el que mejor balanceaba y obtenía el menor resultado de estos 3 grupos de métricas sobre el cual trabajaremos en adelante.

\[ sARIMA(0,1,1)(1,1,2)[2] \]

\[ \Phi_1(B^2)\nabla_1^1 \nabla_2^1 Y_t = \theta_1(B)\Theta(B^2)Z_t\]

Desarrollamos el modelo:
Sea: \[ Y_t=log(X_t) \]  Entonces: \[Y_t = Y_{t-1} + Y_{t-2}(1+\Phi_1) - Y_{t-3}(1+\Phi_1)-\Phi_1Y_{t-4}+\Phi_1Y_{t-5} + Z_t+ \theta_1Z_{t-1}+\Theta_1Z_{t-2}+\theta_1\Theta_1Z_{t-3}+\Theta_2Z_{t-4}+\theta_1\Theta_2Z_{t-5} \]

Coeficientes:

##         ma1        sar1        sma1        sma2 
## -0.87832770 -0.98327076 -0.01056415 -0.84446121

Entonces los parametros son: \[\Phi_1 = -0.98327076,\ \theta_1= -0.87832770, \ \Theta_1= -0.01056415 \ y \ \Theta_2= -0.84446121\]

Intervalos de confianza:

##           2.5 %     97.5 %
## ma1  -1.1409258 -0.6157296
## sar1 -1.0508310 -0.9157105
## sma1 -0.3060396  0.2849113
## sma2 -1.1229379 -0.5659845

Evaluando la significancia de los parámetros encontramos uno no signifcativo, sin embargo, el incluirlo mejora sus propiedades, disminuyendo las métricas del error.

## Series: ts_training 
## ARIMA(0,1,1)(1,1,2)[2] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ma1     sar1     sma1     sma2
##       -0.8783  -0.9833  -0.0106  -0.8445
## s.e.   0.1340   0.0345   0.1508   0.1421
## 
## sigma^2 = 0.1275:  log likelihood = -25.12
## AIC=60.25   AICc=61.4   BIC=70.55
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.01495307 0.4520975 0.3512017 -9.607661 27.87384 0.3793358
##                  ACF1
## Training set 0.094705

Las métricas relacionadas al ajuste del conjunto training son buenas sobre todo, en conjunto, con las demás características que veremos más adelante.

Validación de Supuestos

A continuación, validamos los supuestos de normalidad, homocedasticidad y no correlación sobre los residuales del modelo.

Residuales

No se observan correlaciones significativas.

Estacionariedad

## 
##  Augmented Dickey-Fuller Test
## 
## data:  fit$residuals
## Dickey-Fuller = -3.6908, Lag order = 3, p-value = 0.03304
## alternative hypothesis: stationary

Notemos que P-value < 0.05 por lo tanto rechazamos H0 y tenemos evidencia suficiente para aceptar la hipótesis de que los residuales son estacionarios (A pesar de que más adelante veremos que no se cumple el supuesto de varianza constante).

Normalidad

Comparación de cuantiles de los residuales con los normales

Graficamente, notamos que los datos son en general normales, salvo algunos outliers.

## 
##  Anderson-Darling normality test
## 
## data:  fit$residuals
## A = 0.53422, p-value = 0.165

De la prueba anterior, concluimos que los residuos del modelo son en efecto gaussianos.

Homocedasticidad

Para probar la homocedasticidad, nos centramos en dos pruebas, la Breusch-Pagan Test y la Score Test for Non-Constant Error Variance.

## 
##  studentized Breusch-Pagan test
## 
## data:  Y ~ X
## BP = 3.9309, df = 1, p-value = 0.04741

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 11.19227, Df = 1, p = 0.00082139

Rechazamos ambas pruebas con un nivel de significa del 5%, por lo que tenemos evidencia a favor de que los datos tienen varianza en funcion del tiempo (heterocedasticos). A pesar de no cumplir con este supuesto, decidimos proseguir con su implemetación por dos razones, la primera porque el p-value casi pasa la prueba Breusch-Pagan a un nivel de significancia alpha=5% que es la más potente de entre las dos anteriores y por las métricas de errores que veremos más adelante.

No correlación

Comprobaremos no correlacion en lag’s mas significativos para el numero maximo de lags a considerar, y dado que la serie no presenta un ciclo estacional tan claro; consideramos apropiado usar el criterio de sqrt(n) con n = longitud de datos (training). En nuestro caso resulta apropiado probar con 8 lag’s.

Es decir, para los primeros 8 lags no existe una correlación, por lo cual la hipotesis de independencia de residuales es verosimil.

Sin embargo, es importante recalcar que no podemos considerar rigurosamente que los residuales sean ruido blanco gaussiano pues se falla el supuesto de varianza constante.

Validación de entrenamiento

Ajustamos el modelo propuesto para evaluar las diferencias con los datos de prueba y así obtener una mejor predicción.

Obtenemos las métricas de los errores relacionadas al conjunto de prueba.

  sqrt(mean((forecast(fit, h = 13, lambda = 0)$mean - ts_test)^2))

## [1] 0.3489104

  mean(abs(forecast(fit, h = 13, lambda = 0)$mean - ts_test))

## [1] 0.2760162

Predicción del conjunto de evaluación

Concluimos que el modelo hace un buen ajuste en general, tanto para la parte de entrenamiento como en la predicción del conjunto de evaluación fundamentados en las métricas anteriores y la validación de supuestos, la heterocedasticidad de los residuales provoca un amplio rango en los intervalos de confianza de la predicción, sin embargo, en general es el mejor modelo. Consideremos también que el modelo tuvo el menor AIC de entre todos los modelos probados y todas las transformaciones probadas.

Predicción

En este apartado, hay varios enfoques posibles para hacer el forecasting, en nuestro caso: Reajustamos el modelo obtenido con el conjunto training sobre el total de los datos para poder predecir un margen de 12 años en el futuro a partir de 1925. Esto con el objetivo de que el forecast preserve la tendencia general de la serie pero que a su vez entienda bien la reducción de varianza al final de la serie histórica. Ajustamos modelo elegido con todos los datos históricos.

Se muestran parametros reajustados:

coefficients(modelo_full)

##         ma1        sar1        sma1        sma2 
## -0.89006679 -0.99943606  0.06385804 -0.91250900

Predicción de toda la serie

Predicción en un horizonte de 12 años

Con lo que observamos en la Figure 13 y 14, el modelo asimiló bien la tendencia de los datos originales y mantuvo la varianza observada al final de la serie original. Notemos que el supuesto, usado en el modelo, de que la serie tenga frecuencia 2 ayudó bastante a explicarla a pesar de que en la serie original no había ciclos estacionales significativos. Del mismo modo, el pronóstico resultó apropiado debido a que no hubo un sobreajuste del modelo (aspecto que validamos con las métricas de error en el conjunto training y testing). El ancho de las bandas de confianza es aceptable pese a no cumplir del todo la homocedasticidad. Por lo tanto, nuestro modelo es suficiente y robusto para modelar y precedir correctamente los datos.