Contaminación del aire (+2 años diario)

Serie que muestra la cantidad de partículas contaminantes en el aire \(PM10\), cuya unidad de medida es \({µg}/{m^3}\), se denomina PM10 por sus siglas en inglés “Particulate Matter” a las pequeñas partículas sólidas o líquidas dispersas en la atmósfera, y cuyo diámetro aerodinámico es menor que 10 µm (1 micrómetro corresponde la milésima parte de 1 milímetro).

La serie consta del promedio diario, con muestras tomadas cada hora, desde 2016-08-01 hasta 2018-09-01. Datos obtenidos de SINAICA obtenidos en la estación Benito Juárez (BJU) ubicada en el Valle de México.

#Leer data
library(fpp2)
## Loading required package: ggplot2
## Loading required package: forecast
## Loading required package: fma
## Loading required package: expsmooth
library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(ggplot2)
library(gridExtra)
data <- read.csv('/home/israel/Documentos/R/Econometria EARF/SeriePM10.csv')
PM10=data$Valor; #data$Valor=data[,2]
inds <- seq(as.Date("2016-08-01"), as.Date("2018-09-01"), by = "day")
tsPM10 <- ts(PM10,
             start = c(2016, as.numeric(format(inds[1], "%j"))),
             frequency = 365)
plot(tsPM10, col='blue',main='PM10',xlab="",ylab='µg/m^3')

Como se puede apreciar en el grafico se puede intuir que existe autocorrelación y no muestra alguna clase de pendiente en su tendencia, pero si muestra una frecuencia o periodo.

Descomposición de series de tiempo

Método de medias móviles

G1 <- autoplot(tsPM10, series="PM10") +
  autolayer(ma(tsPM10,3,centre=F), series="3-MA") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")

G2 <- autoplot(tsPM10, series="PM10") +
  autolayer(ma(tsPM10,5,centre=F), series="5-MA") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")

G3 <- autoplot(tsPM10, series="PM10") +
  autolayer(ma(tsPM10,7,centre=F), series="7-MA") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")

G4 <- autoplot(tsPM10, series="PM10") +
  autolayer(ma(tsPM10,9,centre=F), series="9-MA") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")

grid.arrange(G1,G2,G3,G4)
## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 6 rows containing missing values (geom_path).
## Warning: Removed 8 rows containing missing values (geom_path).

Como se puede apreciar en los procesos autoregresivos 3, 5, 7 y 9 no suaviza la serie.

startdate <- "2016-08-01" 
enddate <- "2018-09-01" 
t= as.numeric(as.Date(startdate)):as.numeric(as.Date(enddate))
p=periodogram(PM10)

dd = data.frame(freq=p$freq, spec=p$spec)
order = dd[order(-dd$spec),]
top3 = head(order, 3)
top3
##          freq      spec
## 2 0.002604167 66812.837
## 3 0.003906250  8078.082
## 7 0.009114583  7334.626

Para detectar una estacionalidad, detectando la fuerza de cada frecuencia posible, y se observa que la mas grande esta en los 0.0026 Hz.

time = 1/top3$f
time
## [1] 384.0000 256.0000 109.7143

Como la serie es diaria, el tiempo en el que se muestra un ciclo es a los 384 días, por falta de datos no podemos completar la suficiente información para que se muestre que el ciclo es cada 365 dias.

autoplot(tsPM10, series="PM10") +
  autolayer(ma(tsPM10,384,centre=F), series="384-MA") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")
## Warning: Removed 383 rows containing missing values (geom_path).

Al aplicar una media móvil de 384 días, se ve clara una serie suavizada, pero se pierde la mitad de la información.

Descomposición clásica

fit <- decompose(tsPM10, type='additive')
autoplot(fit)

autoplot(tsPM10, series="tsPM10") +
  autolayer(seasadj(fit), series="Seasonally adj. data") +
  xlab("Year") + ylab("µg/m^3") +
  ggtitle("Particulate Matter 10 µm: Valle de México")

En la descomposición clásica no ayuda para simplificar la serie, por lo que se utilizara un modelo automático.

autoplot(tsPM10) + geom_smooth(method="auto", se=FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

En este modelo se muestra un polinomio que no parece ser el más adecuado. Como no muestra una tendencia que sea suave y conociendo que la información es anual, Se proponen modelos trigonométricos y se verifica si hay una tendencia lineal.

Modelo Trigonométrico

ssp <- spectrum(PM10)

per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]

rg <- diff(range(PM10))
plot(PM10~t,ylim=c(min(PM10)-0.1*rg,max(PM10)+0.1*rg),col=8)

modlin <- lm(PM10~t)
summary(modlin)
## 
## Call:
## lm(formula = PM10 ~ t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.238 -12.358  -2.032  10.489  55.903 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.3045397 45.9229465   1.204    0.229
## t           -0.0008375  0.0026399  -0.317    0.751
## 
## Residual standard error: 16.03 on 760 degrees of freedom
## Multiple R-squared:  0.0001324,  Adjusted R-squared:  -0.001183 
## F-statistic: 0.1006 on 1 and 760 DF,  p-value: 0.7512
lines(fitted(modlin)~t,col=2,lty=3,lwd=2.5)

reslm <- lm(PM10 ~ sin(2*pi/per*t)+cos(2*pi/per*t))
summary(reslm)
## 
## Call:
## lm(formula = PM10 ~ sin(2 * pi/per * t) + cos(2 * pi/per * t))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.773  -8.898  -0.628   8.187  49.827 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          40.6326     0.4704  86.377  < 2e-16 ***
## sin(2 * pi/per * t) -13.0423     0.6674 -19.541  < 2e-16 ***
## cos(2 * pi/per * t)   2.8373     0.6631   4.279 2.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.98 on 759 degrees of freedom
## Multiple R-squared:  0.3448, Adjusted R-squared:  0.3431 
## F-statistic: 199.7 on 2 and 759 DF,  p-value: < 2.2e-16
lines(fitted(reslm)~t,col=4,lty=1,lwd=2.5)

reslm2 <- lm(PM10 ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
summary(reslm2)
## 
## Call:
## lm(formula = PM10 ~ sin(2 * pi/per * t) + cos(2 * pi/per * t) + 
##     sin(4 * pi/per * t) + cos(4 * pi/per * t))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.243  -8.464  -0.636   7.378  46.683 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          40.6216     0.4629  87.755  < 2e-16 ***
## sin(2 * pi/per * t) -13.0634     0.6568 -19.889  < 2e-16 ***
## cos(2 * pi/per * t)   2.8440     0.6525   4.359 1.49e-05 ***
## sin(4 * pi/per * t)   3.3311     0.6537   5.096 4.39e-07 ***
## cos(4 * pi/per * t)  -0.6200     0.6556  -0.946    0.345    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.78 on 757 degrees of freedom
## Multiple R-squared:  0.3673, Adjusted R-squared:  0.364 
## F-statistic: 109.9 on 4 and 757 DF,  p-value: < 2.2e-16
lines(fitted(reslm2)~t,col=3,lty=2,lwd=2.5) 

reslm3 <- lm(PM10 ~ 
               sin(2*pi/365*t)+cos(2*pi/365*t)+
               sin(2*pi/7*t)+cos(2*pi/7*t)+
               cos(4*pi/365*t))
summary(reslm3)
## 
## Call:
## lm(formula = PM10 ~ sin(2 * pi/365 * t) + cos(2 * pi/365 * t) + 
##     sin(2 * pi/7 * t) + cos(2 * pi/7 * t) + cos(4 * pi/365 * 
##     t))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.940  -7.880  -0.342   7.052  43.741 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          41.3544     0.4383  94.345  < 2e-16 ***
## sin(2 * pi/365 * t)  10.4044     0.6154  16.907  < 2e-16 ***
## cos(2 * pi/365 * t)   8.5324     0.6243  13.668  < 2e-16 ***
## sin(2 * pi/7 * t)    -1.9298     0.6185  -3.120  0.00188 ** 
## cos(2 * pi/7 * t)     3.1904     0.6190   5.154 3.26e-07 ***
## cos(4 * pi/365 * t)   4.3793     0.6267   6.987 6.15e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.08 on 756 degrees of freedom
## Multiple R-squared:  0.4354, Adjusted R-squared:  0.4317 
## F-statistic: 116.6 on 5 and 756 DF,  p-value: < 2.2e-16
lines(fitted(reslm3)~t,col=1,lty=1,lwd=2)

En el grafico se muestran: 1. Serie original en puntos grises 1 2. Linea Tendencia lineal de color Rojo (se observa que tiene una pendiente negativa pero no es significativa) 2 3. Primer aproximación trigonométrica donde el tiempo depende de Senos y Cosenos de color Azul 3 4. Segunda aproximación Trigonométrica agregando una segunda armonía de color Verde 4 5. Tercer aproximación Trigonométrica suponiendo 2 términos periódicos, anual y semanal, de color Negro. Que es la función a utilizar que es la que nos muestra una mejor tendencia. 5

Análisis de los residuos y correlograma

par(mfrow = c(2, 2))
plot(y=rstandard(reslm3),x=as.vector(time(PM10)),xlab='Tiempo',
     ylab='Residuales estandarizados',type='o')

hist(rstandard(reslm3),xlab='Residuales estandarizados')

qqnorm(rstandard(reslm3));qqline(rstandard(reslm3), col='red')

acf(rstandard(reslm3), lag.max=31)

En la primer grafica se muestran los residuos, lo que nos sugiere que aun falta limpiar la serie para que no se muestre alguna dependencia. En el histograma se puede apreciar que los residuos están pareciéndose a una normal estándar junto con la grafica Q-Q, aunque en los extremos aun no se ajustan a una distribución normal y por último en el correlograma se puede observar que hay varias dependencias autoregresivas muy marcadas en los primeros 14 ordenes.