Proyecto Series de Tiempo

Proyecto 2: Series de Tiempo

Dirección de Monitoreo Atmosférico

Introducción

Entre los procesos que originan la contaminación del aire en nuestra ciudad se encuentran las partículas suspendidas, que son cualquier tipo de material sólido o líquido que se encuentra en suspensión en el aire. En la Ciudad de México una fracción importante se forma de reacciones químicas en la atmósfera contaminada aunque entre otras fuentes de emisión de este contaminante se encuentran las industrias, los incendios, las emisiones de camiones y automóviles. Su tamaño varía y va desde unas cuantas millonésimas de milímetros hasta llegar a alcanzar el tamaño de un grano de arena.

Para el presente proyecto se tendrán en cuenta solo dos de las partículas suspendidas, PM10 y PM2.5, la diferencia entre ellas radica en que el diámetro de las PM10 es inferior a 0.01 milímetros y las PM2.5 tienen un diámetro inferior a 0.0025 milímetros, son tan finas que no es posible verlas a simple vista.

Las partículas suspendidas representan el principal problema de salud pública, ya que agravan enfermedades cardiovasculares y respiratorias como el asma. Además de que la exposición crónica a altas concentraciones puede provocar un incremento en el riesgo de morbilidad y mortalidad.

Sabemos que la población esta creciendo cada vez más y así como más personas estan expuestas a los riesgos de salud que estas partículas representan, tambíen participan en la contaminación de las mismas, por lo que es necesario realizar un análisis a través del tiempo de la cantidad de partículas suspendidas en la ciudad para tener bases con las que se pueda trabajar una propuesta a partir de un pronóstico dado.

El objetivo será entender el comportamiento de a concentracion de las particulas en los proximos años, y usar nuestras predicciones de la serie de tiempo para saber si la calidad del aire para el Valle de México empeorará o mejorará.

Descripción de la base de datos

Los datos con los que se trabajarán son estadísticas registradas para la Ciudad de México y área metropolitana desde enero del 2004 a octubre del 2019, debido a que apartir del 2004 se consideraron las mediciones de las particulas PM2.5.

Nuestra base de datos cuenta con 5 variables: una referida a la fecha, la segunda al nombre del identificador, si este registra partiulas PM10 o PM2.5, el valor que registró el identificador de particulas, así como la escala de medición, de las cuales las primeras tres son cualitativas y las ultimas dos numéricas. El número de registros presentes en esta base es de 143,542.

Cabe recalcar que de las 5 variables solo utilizamos 3 ya que al colapsar las mediciones diarias a mensuales no nos importa el lugar exacto de la medición de cada día. Tampoco consideramos para el análisis la columna unit porque esta tiene la misma medición para todas las particulas, es decir de microgramos por metro cúbico.

Entonces las variables usadas son las siguientes:

• date: el día en que se tomó la medición.

• id_parameter: el tipo de partícula que se midió, ya sea PM10 o PM2.5.

• value: el promedio de partículas medidas en ese mes.

Análisis exploratorio de datos

Gráficas de dispersión y relación de los datos, mostrando tambien la correlación

Se observa que las partículas PM2.5 no presentan una tendencia a través del tiempo ni alguna concentración, es decir, se encuentran dispersos por toda la gráfica, también podemos diferenciar el caso del día 16 de mayo donde la concentración de las partículas PM2.5 fue muy alta. Hablando de la dispersión de los datos, es similar el caso de las partículas PM10 por lo que no se incluye su gráfico. En la gráfica de Relación entre partículas se observa una tendencia lineal creciente por lo que se considerará una regresión lineal y también se calculará la correlación entre ellas.

[1] 0.9121002

Vemos que la correlación entre las partículas es muy alta; ajustando un modelo de regresión lineal tenemos que por cada unidad de PM2.5 que se aumente en la CDMX y área metropolitana, las particulas PM10 aumentarán en 2.027 unidades y la variable PM2.5 describe significativamente a la variable PM10. Vale la pena analizar las dos partículas ya que varían en sus fechas y lugar de medición por lo que los resultados podrán ser diferentes.

Los histogramas y densidades son:

Es claro que los datos se concentran en menor cantidad para valores altos ya que el gobierno implementa medidas para controlar el aumento de las partículas PM2.5 y PM10 con programas como “No circula” cuando dichas partículas incrementan.

Veremos en mapas el cambio de los valores de las partículas en los años 2005 y 2019; el color azul indica que no hubo mediciones durante ese año en dicho municipio, los círculos indican a través de su tamaño y color, que tan elevado es el valor de dicha partícula, entre más grande y más rojo significa que es mayor el nivel de contaminación por culpa de la partícula.

El siguiente mapa es de las partículas PM10 en el año 2005

El siguiente mapa es de las partículas PM10 en el año 2019

El siguiente mapa es de las partículas PM2.5 en el año 2005

El siguiente mapa es de las partículas PM2.5 en el año 2019

Para ambos casos se observa que las partículas contaminantes incrementaron en todas las zonas.

Modelo para particulas pm10.

Primeramente veamos como es la serie de tiempo para esta partícula a lo largo del tiempo.

A simple vista se puede notar que esta serie no será completamente estacionaria dado que, aunque la media si parece ser constante, la varianza, al variar demasiado puede afectar este aspecto de las hipótesis de una serie estacionaria; para estar completamente seguros de esto planteremos tres pruebas de hipótesis con las cuales podremos concluir con evidencia estadística la observación anterior:

• adf.test(): Una prueba de Dickey-Fuller aumentada (ADF) es una prueba de raíz unitaria para una muestra de una serie de tiempo. La estadística Dickey-Fuller Aumentada (ADF), utilizada en la prueba, es un número negativo. Cuanto más negativo es, más fuerte es el rechazo de la hipótesis nula de que existe una raíz unitaria para un cierto nivel de confianza. Las hipótesis son:
  H0: La serie no es estacionaria vs H1: La serie es estacionaria

• kpss.test(): La hipótesis nula para la prueba KPSS es que los datos son estacionarios. Para esta prueba, NO queremos rechazar la hipótesis nula. En otras palabras, queremos que el p-valor obtenido sea mayor que 0.05, no menor que 0.05. Las hipotesis son: H0: Los datos son estacionarios vs H1: Los datos no son estacionarios

    Augmented Dickey-Fuller Test

data:  ts_pm10
Dickey-Fuller = -9.3557, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

    KPSS Test for Level Stationarity

data:  ts_pm10
KPSS Level = 0.53592, Truncation lag parameter = 4, p-value = 0.03358

De lo anterior:

Para adf.test se tiene: p-value = 0.01 < 0.05, entonces podemos aceptar que sea estacionaria. Sin embargo para kpss.test se tiene: pvalue = 0.03358 < 0.05, entonces la serie es no estacionaria. Dado que una prueba falló la serie no es estacionaria.

Revisemos las funciones de autocorrelación de las observaciones para dar una idea de las dependencias entre las observaciones:

De la gráfica anterior en cuanto a las correlaciones totales se aprecia una gran estacionalidad en la serie, mientra que la parcial se comporta mucho mejor.

Para remediar ésto se suavizara la serie por medio de medias moviles, con el fin de disminuir la varianza de las observaciones.

Entonces la gráfica anterior se nota mucho más manejable a comparación de la original, así que se trabajará con esta.

Analizando los correlogramas se observa que los problemas se mantienen, entonces se aplicará una descompocisión de la serie en sus componentes, y se trabajará con una serie libre del factor tendencial y estacionario, es decir, se trabajará con la parte aleatoria de está.

De ésto, la descompocisión queda como:

Asimismo la serie de tiempo ts10 es aquella a la cual se le omitieron estas componenetes

Chequemos las autocorrelaciones de esta serie:

Entonces las autocorrelaciones totales se comportan de mejor manera que la serie original y aparentente ya se comporta de una manera más estacional.

Ahora revisaremos si esta serie es estacionaria con las pruebas de hipótesis vistas antes.

    Augmented Dickey-Fuller Test

data:  na.omit(ts10)
Dickey-Fuller = -5.4892, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

    KPSS Test for Level Stationarity

data:  na.omit(ts10)
KPSS Level = 0.034649, Truncation lag parameter = 4, p-value = 0.1

Para esta serie se aprecia que bajo la prueba Dickey-Fuller se tiene: p-value = 0.01 < 0.05, considerando así una serie estacionaria; de igual modo, para la prueba KPSS se obtuvo: p-value = 0.1 > 0.05, no rechazando la hipotesis nula, asumiendo de este modo una serie estacionaria.

Dado que ya se obtuvo una serie que es estacionaria se puede obtener el modelo que ajuste de buena manera las observaciones, cabe recordar que es un modelo para la parte residual de la serie original suavizada.

Se hará uso de la función auto.arima cargada en Rstudio.

Series: ts10 
ARIMA(2,0,1)(1,0,0)[12] with zero mean 

Coefficients:
         ar1     ar2     ma1     sar1
      0.0023  0.5230  0.9613  -0.2145
s.e.  0.0853  0.0819  0.0407   0.0811

sigma^2 estimated as 1.43:  log likelihood=-275.76
AIC=561.51   AICc=561.87   BIC=577.28

De modo que el mejor modelo para esta serie consta, de un modelo autorregresivo de orden 2, tomando las dos observaciones anteriores a la que se encuentra en tiempo t, de un parámetro en medias moviles tomando error anterior a la observación en tiempo t, además consta de un parámetro en la parte estacionaria autoregresivo, es decir, se considera para el modelo la observación regisrada en el tiempo t-12.

Ahora, con el modelo ajustado, comparemos los datos reales con los obtenidos por el modelo:

Se observa que los residuales prensentan un comportamiento, en media, alrededor de cero mientras que la varianza no varía mucho entre ellos.

Predicciones

Bajo este modelo obtengamos las predicciones a 18 meses.

Notamos que esta predicción converge a la media relativamente rápido, por lo que, este modelo solo sería prático para predicciones a corto plazo.

Aun así busquemos regresar al modelo suavizado original:

A partir de la linea vertical (hacia la derecha) se encuentran los valores de la predicción para el modelo a pártir de mayo de 2019 hasta principios de 2020.

Se calcula una tasa de error relativo para conocer si nuestro modelo es poco equívoco.

[1] 0.01540677

Se observa que el modelo se equivocó un 1.5% al predecir los valores con los que ya contábamos.

Modelo para particulas pm2.5

Veamos como son los datos para esta serie de tiempo:

Al igual que la serie anterior, evidentemente hay muchos problemas en cuanto a la varianza de la serie.

Para estas particulas se hizo un procedimiento análogo al modelo anterior, llegando a una serie que consta solo del componenente aleatorio (ts25) de la serie y que es estacionaria, del cual veremos las gráficas de autocorrelación respectivas.

Apreciamos que la autocorrelación total se debe corregir un poco mientras que la parcial no presenta gran problema.

Para probar que esta serie es estacionaria haremos las mismas pruebas de hipótesis que anteriormente se hicieron para ver si efectivamente es un serie estacionaria.

    Augmented Dickey-Fuller Test

data:  ts25
Dickey-Fuller = -7.0851, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

    KPSS Test for Level Stationarity

data:  ts25
KPSS Level = 0.052284, Truncation lag parameter = 4, p-value = 0.1

Para ambas pruebas se obtuvo que no se puede rechazar el hecho de considerar la serie como estacionaria, por lo que el primer objetivo, obtener una serie estacionaria, ya se obtuvo; entonces se ajustará un modelo usando autoarima de Rstudio

Series: ts25 
ARIMA(0,0,4)(2,0,0)[12] with zero mean 

Coefficients:
         ma1     ma2      ma3      ma4     sar1     sar2
      1.2477  0.6928  -0.0118  -0.4031  -0.1509  -0.1536
s.e.  0.0927  0.1193   0.1427   0.1107   0.0923   0.0879

sigma^2 estimated as 0.3186:  log likelihood=-145.92
AIC=305.84   AICc=306.52   BIC=327.91

El mejor ajuste de la parte aleatoria se ajusta son un SARIMA donde la parte de medias moviles se ajusta los 4 errores anteriores a la observacion en t, mientras que la parte de la estacionalidad consta de dos parámetros pero autorregresivos con saltos cada 12 periodos, es decir, t-12 y t-24.

Comparemos el modelo suavizado con el modelo ajustado:

De las anteriores tres gráficas: las primeras dos intuyen a pensar que los valores de la serie que se suavizó y los ajustados concuerdan en gran medida, mientras que el último gráfico son los residuales del modelo.

Predicciones

Bajo este modelo obtengamos las predicciones a 18 meses.

Las predicciones obtenidas para este modelo se aprecian que convergen a la media, por lo que a largo plazo este modelo podría no ser muy efectivo si se quieren proyecciones a periodos de tiempo grandes, este caso sucede al menos con los errores del modelo.

Ahora, estas predicciones son para el modelo respecto a los errores ajustados, regresemos al modelo, sumando la tendencia y estacionalidad respectiva para poder ver un aproximado de predicciones para la serie suavizada original.

Los puntos que se encuntran a la derecha de la línea vertical son las predicciones bajo este modelo para las mediciones de partículas.

Se calcula una tasa de error relativo para conocer si nuestro modelo es poco equívoco.

[1] 0.06121523

Se observa que el modelo se equivocó un 6.1% al predecir los valores con los que ya contábamos.

Podemos concluir que nuestros modelos son adecuados para predecir los futuros valores de las partículas PM2.5 y PM10 en el Valle de México con pequeñas tasas de error, un problema que dado que las predicciones convergen a la media relativamente rápido no es conveniente predecir a largo plazo (más de año y medio), por lo que los errores aumentarían.