MODELOS ARIMA

MODELOS ARIMA

1. Breve descripciòn de la serie elegida, asi como su frecuencia, unidad de medida, fuente, etc.

library(forecast)
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(urca)
library(ggplot2)

library(foreign)
prodpesquera= read.csv("~/SEPTIMO SEMESTRE/SERIES DE TIEMPO/1ER PARCIAL/prodpesquera.csv")
class(prodpesquera)

## [1] "data.frame"

• Información de la base de datos:
  • Ruta temática:
    • Sector alimentario> Producción> Captura pesquera> Volumen de la captura pesquera en peso vivo por litoral y entidad federativa> Total
  • Periodicidad:
    • Mensual
  • Unidad de medida:
    • Toneladas
  • Fuente:
    • SAGARPA. Comisión Nacional de Acuacultura y Pesca.
  • Fecha inicial:
    • 2015/12
  • Fecha final:
    • 2002/01
  • Última actualización:
    • 2016/06/27 DATOS DE PRODUCCIÓN PESQUERA

Enlace en línea

*Fijar la base de datos.

attach(prodpesquera)

*Abrir la base de datos.

View(prodpesquera)

Se declara como “Serie de Tiempo”

pp<-ts(prodpesquera$Dato, frequency = 12, start = c(2002, 1))
class(pp)

## [1] "ts"

2. Graficar los datos, analizar patrones y observaciones, atìpicas. Apoyen su analisis con una descomposicion clàsica.

autoplot(pp, colors="orange", main= "PRODUCCIÓN PESQUERA EN MÉXICO", ylab="toneladas", xlab= "años")

En la gráfica “PRODUCCIÓN PESQUERA EN MÉXICO”, Se pueden analizar los componentes de la serie de tiempo. Para este caso, la serie presenta estacionalidad y tendencia-ciclo, para poder apreciar mejor este hecho, se analizara la grafica por sus componentes.

En la gráfica de arriba se muestran los datos originales, comprobando la estacionalidad de serie de tiempo y los ciclos que muestra. a continuación se muestran los componentes de la serie:

Descomposición de la serie de tiempo

Fit<-decompose(pp)

library(ggplot2)
plot(Fit, col= "pink")

• Tendencia-ciclo:
  • Permite apreciar el comportamineto de la producción pesquera en el pais a través del tiempo. su tendencia es lineal aunque actualmente la producción pesquera esta disminuyendo (segun las cooperativas del pais) y aparenta ciclos ya que varian los años en que aumenta y en que cae la producción precipitadamente.
  • Esta variación en la producción se debe a muchos aspectos, el mas representativo es el factor “lluvia” ya que si en algún año se presentan sequías, es muy probable que la produción pesquera sea baja.
  • Al no haber lluvias se baja considerablemente el nivel de agua en la zona de esteros y pampas, situación por la que la producción disminuye drásticamente, pegándole a los pescadores que ven impactada su economía.
  • A LA BAJA LA PRODUCCIÓN PESQUERA DEL PAÍS

Gráfica de estacionalidad #1

ggseasonplot(pp, year.labels = TRUE, year.labels.left = TRUE)+ ylab("Toneladas")+ xlab("Años") + ggtitle("PRODUCCIÓN PESQUERA EN MÉXICO")

Se pueden apreciar mejor los meses que tienen estacionalidad, por ejemplo; septiembre y febrero.

Gráfica de estacionalidad #2

ggseasonplot(pp, Main= ("PRODUCCIÓN PESQUERA EN MÉXICO"), polar="TRUE") 

Con esta gráfica se puede observar mejor la estacionalidad, mostrando claramente que en febrero, abril y septiembre se tiene una menor producción pesquera en el país.

3. Si es necesario, utilizar transformaciòn Box-Cox/logaritmos para estabilizar la varinaza.

Transformaciòn Box-Cov

BoxCox.ar(pp)

El intervalo a 95% para \(\lambda\), contiene un valor de \(\lambda=0\) muy cerca de su centro y sugiere una transformación logarítmica.

pplog <- log(pp)
plot(pplog, main="PROD.PESQUERA EN MÉXICO", col="pink", ylab="toneladas", xlab="años", lwd=2, type="o", pch=1)

Para estabilizar la varianza se aplican logarítmos (comprime los datos). sin embambargo la serie sigue teniendo problemas con la media. asi que le quitaremos la tendencia con diferencias.

4.En caso de estacionalidad, aplicar diferencias estacionales.

serie Logarítmica

pplog <- log(pp)
autoplot(pplog, main="PROD.PESQUERA EN MÉXICO", col="pink", ylab="toneladas", xlab="años", lwd=2, type="o", pch=1)

Primera diferencia

difpp <- diff(pplog)
plot(diff(log(pp)), main='Primera Diferencia de Prod. Pesquera en México')

Segunda Diferencia

dif2pp <- diff(pp,difference=2)
plot(diff(pp,difference=2), main='Segúnda Diferencia de Prod. Pesquera en México')

library(gridExtra)
p1<- autoplot (log(pp)) + xlab("años") + ylab("toneladas") + ggtitle("Log. de PROD.PESQUERA EN MÉXICO")
p2<- autoplot (diff(log(pp))) + xlab("años") + ylab("toneladas") + ggtitle("Primera Diferencia de PROD.PESQUERA EN MÉXICO")
p3<- autoplot (diff(pp, difference=2)) + xlab("años") + ylab("toneladas") + ggtitle("Segunda Diferencia de PROD.PESQUERA EN MÉXICO")

grid.arrange(p1, p2, p3)

Ya con las diferencias la media se vuelve estable, alrededor de 0.

*Una diferencia de logaritmos puede ser vista tambien como una tasa de crecimeinto.

Tasa de crecimeinto de la Producción pesquera en México.

tasa<-na.omit((pp-zlag (pp))/ zlag(pp))
plot(x=diff(log(pp))[-1], y=tasa[-1])

La gráfica nos muestra una correlación positiva, ya que la recta es muy cercana a 45 grados, eso quiere decir estan saliendo casi los mismos datos. Y eso lo podemos checar con una correlación.

cor(diff(log(pp))[-1], tasa[-1])

## [1] 0.9893194

La correlación es de 0.98; muy cercana a 1. El coeficiente de correlación comprueba lo ya visto en la gráfica; que la diferencia de los logarímos es practicamente una tasa de crecimiento común.

5.Usar prueba Dickey-Fuller para evaluar el orden de integración de la serie. diferenciar hasta que la serie sea estacionaria.

summary(ur.df(dif2pp, type = "none", selectlags = "AIC"))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -112442  -22026    3057   21515   97715 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.97691    0.12329  -16.04  < 2e-16 ***
## z.diff.lag  0.33578    0.07175    4.68 5.74e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32690 on 174 degrees of freedom
## Multiple R-squared:  0.7675, Adjusted R-squared:  0.7648 
## F-statistic: 287.2 on 2 and 174 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -16.0346 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Para que la serie sea estacionaria: phi tiene que ser menor a 1.

H0=No es estacionaria. H1=Ya es estacionaria.

Al tomar la segunada diferencia, se muestra que la serie ya es estacinaria, ya que el valor estadistico es de -16.0346, es mayor a los valores críticos de Tau.

Una caminata aleatoria no es estacionaria. Un ruido blanco si es estacionario.

6.Usar herramientas ACF, PACF, EACF y criterios de Akaike/Bayes para construir propuestas de modelos.

ACF (MA)

Se van a utilizar los datos ya transformados y estacionarios.

ggAcf(dif2pp)

Se propone un MA(1)…

PACF (AR)

ggPacf(dif2pp) 

Se propone un AR(2)…

EACF (ARMA)

eacf(dif2pp)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o x x o o  x  o  o 
## 1 x o o o o o o x o o o  x  o  o 
## 2 x x x o o o o o o o o  x  o  o 
## 3 x o x o o o o o o o o  x  o  o 
## 4 x x x o x o o o o o o  x  o  o 
## 5 x o x o x o o o o o o  x  o  o 
## 6 x o x o o o o o o o o  x  o  o 
## 7 o x x x o x x x o o o  x  o  o

Esta gráfica no es clara para esta serie, así que no ayuda mucho…

ggtsdisplay(dif2pp)

Se observan las opciones de MA(1) y AR(2).

Criterio de Akaike-Bayes.

res <- armasubsets(y=dif2pp,nar=10,nma=10,y.name='test',ar.method='ols') 

plot(res)

La gráfica muestra que puede ser un AR(1) o un MA(1).

Propuesta #1 AR(2)

propuesta1<-arima(dif2pp,order = c(2,0,0)) 
propuesta1

## 
## Call:
## arima(x = dif2pp, order = c(2, 0, 0))
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.6380  -0.3321    33.5283
## s.e.   0.0708   0.0705  1234.6905
## 
## sigma^2 estimated as 1.045e+09:  log likelihood = -2101.14,  aic = 4208.28

El criterio de akaike es muy alto.

Propuesta #2 MA(1)

propuesta2<-arima(dif2pp,order = c(0,0,1)) 
propuesta2

## 
## Call:
## arima(x = dif2pp, order = c(0, 0, 1))
## 
## Coefficients:
##          ma1  intercept
##       -1.000    -4.2035
## s.e.   0.014    37.6272
## 
## sigma^2 estimated as 674354740:  log likelihood = -2064.47,  aic = 4132.93

El criterio de akaike nos dice que mejoró el modelo al disminuir el criterio de akaike.

ANALIZAR LOS RESIDUOS…

Residuos de propuesta #1 ar(2)

ggAcf(residuals(arima((dif2pp),order=c(2,0,0))))

Los residuos muestran que existe estacionalidad, por tal motivo un modelo sólo No estacionario no es suficiente para demostrar el comportamiento de la serie que se pretende modelar.

Residuos de propuesta #2 MA(1)

ggAcf(residuals(arima((dif2pp),order=c(0,0,1))))

Las gráficas de los residuos para el caso de ambas propuestas muestran que la serie tiene estacionalidad por lo que…

checkresiduals(propuesta1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0) with non-zero mean
## Q* = 159.33, df = 21, p-value < 2.2e-16
## 
## Model df: 3.   Total lags used: 24

checkresiduals(propuesta2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,1) with non-zero mean
## Q* = 116.02, df = 22, p-value = 9.104e-15
## 
## Model df: 2.   Total lags used: 24

Al parecer el modelo es un ARIMA estacional. Por lo que el modelado se modelará otra vez..

MODELO ARIMA ESTACIONAL

Un modelo ARIMA ESTACIONAL se forma:

\((p,d,q)\) No Estacional. \((P,D,Q)_s\) Estacional.

La diferencia estacional del periodo \(s\) para la serie \({Y_t}\) se denota \(\nabla_s\) \(Y_t=Y_t-Y_{t-s}\)

Se obtienen los logaritmos

Existe un pequeño aumento en la varianza, así que se tomaran logarítmos para estabilizar la varianza.

pplog <- log(pp)
autoplot(pplog)+ylab('log(Toneladas)')+xlab('años')+ggtitle('Prod. Pesquera en México')

Primer Diferencia Típica

ggtsdisplay(difpp, main= ("Primer Dif Tipica del Log de Prod. Pesquera"))

Existe un pequeño aumento en la varianza, así que se tomaran logarítmos para estabilizar la varianza.

Segunda Diferencia Típica

ggtsdisplay(dif2pp, main= ("Segunda Dif Tipica del Log de Prod. Pesquera"))

Existe un alto componente estacional por lo que se aplica diferencia estacional.

Primer Diferencia Estacional

ggtsdisplay(diff(pplog, 12), main= ("Primer Dif estacional del Log de Prod. Pesquera"))

El pico significativo en el rezago 1 del PACF sugiere un componente AR(1) no estacional, y el pico significativo en el rezago 12 sugiere un componente estacional AR(1) En los gráficos de la diferencia estacionaria, hay patrones parecidos en AFC y PACF en los rezagos 1, 12, y vecinos.

Propuesta ARIMA(1,0,1)(1,1,1) #1

Esto sugiere iniciar con un modelo ARIMA (1,0,1)(1,1,1)

fit1<- arima(pplog, order=c(1,0,1), seasonal=c(1,1,1))
fit1

## 
## Call:
## arima(x = pplog, order = c(1, 0, 1), seasonal = c(1, 1, 1))
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.7428  -0.2946  0.0684  -0.8273
## s.e.  0.1000   0.1504  0.1136   0.0970
## 
## sigma^2 estimated as 0.0175:  log likelihood = 95.05,  aic = -182.09

Propuesta ARIMA(2,0,1)(1,1,2) #2

Ahora se propondrá con un modelo ARIMA (2,0,1)(1,1,2)

fit2<- arima(pplog, order=c(2,0,1), seasonal=c(1,1,2))
fit2

## 
## Call:
## arima(x = pplog, order = c(2, 0, 1), seasonal = c(1, 1, 2))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sma1     sma2
##       1.2936  -0.3357  -0.8239  -0.5031  -0.2752  -0.4394
## s.e.  0.1642   0.1265   0.1283   0.9619   0.9483   0.7401
## 
## sigma^2 estimated as 0.01734:  log likelihood = 96.01,  aic = -180.01

Se mostraran los residuos para validar confirmar la validación del modelo.

7.Como parte del diagnostico del modelo analizar los residuos y aplicar la prueba Ljung-Box.

Residuos Propuesta Arima #1

fit1<- arima(pplog, order=c(1,0,1), seasonal=c(1,1,1))
fit1

## 
## Call:
## arima(x = pplog, order = c(1, 0, 1), seasonal = c(1, 1, 1))
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.7428  -0.2946  0.0684  -0.8273
## s.e.  0.1000   0.1504  0.1136   0.0970
## 
## sigma^2 estimated as 0.0175:  log likelihood = 95.05,  aic = -182.09

checkresiduals(fit1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(1,1,1)[12]
## Q* = 11.164, df = 20, p-value = 0.9419
## 
## Model df: 4.   Total lags used: 24

Con el valor “P-VALUE” se comprueba la hipotesis de los residuos agrupados. H0:(+0.05) No correlación en los residuos grupales.

H1:(-0.05) Hay correlación serial en los residuos.

En la primer propuesta se obtiene un p-value=0.94, lo que indíca que no hay correlación serial.

Residuos Propuesta Arima #2

fit2<- arima(pplog, order=c(2,0,1), seasonal=c(1,1,2))
fit2

## 
## Call:
## arima(x = pplog, order = c(2, 0, 1), seasonal = c(1, 1, 2))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sma1     sma2
##       1.2936  -0.3357  -0.8239  -0.5031  -0.2752  -0.4394
## s.e.  0.1642   0.1265   0.1283   0.9619   0.9483   0.7401
## 
## sigma^2 estimated as 0.01734:  log likelihood = 96.01,  aic = -180.01

checkresiduals(fit2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(1,1,2)[12]
## Q* = 11.183, df = 18, p-value = 0.8864
## 
## Model df: 6.   Total lags used: 24

H0:(+0.05) No correlación en los residuos grupales.

H1:(-0.05) Hay correlación serial en los residuos.

En la primer propuesta se obtiene un p-value=0.88, lo que indíca que no hay correlación serial. Tambien se acepta la Hipotesis nula.

En conclusión es evidente que los residuos individuales y grupales quedaron mucho mejor con la propuesta ARIMA #1, ya que el criterio de akaike es menor, en la prueba de Ljung-Box es de mayor valor y lo mejor es que los residuos individuales y grupales ya no tienen correlación.

8.Usar la función auto.arima() y compare sus resultados con el inciso anterior

propuesta.arima<-auto.arima(pp, stepwise = FALSE, approximation = FALSE)
propuesta.arima

## Series: pp 
## ARIMA(2,0,1)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       1.3213  -0.3539  -0.8575  -0.8010
## s.e.  0.1329   0.1071   0.0961   0.0746
## 
## sigma^2 estimated as 333244907:  log likelihood=-1890.8
## AIC=3791.6   AICc=3791.97   BIC=3807.22

A continuación se le aplicará la prueba Ljung-Box y se graficaran los residuos de la propuesta arrojada por “Auto.arima”

Propuesta.arima <- arima(pp, order=c(2,0,1), seasonal=c(0,1,1))
Propuesta.arima

## 
## Call:
## arima(x = pp, order = c(2, 0, 1), seasonal = c(0, 1, 1))
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       1.3213  -0.3539  -0.8575  -0.8010
## s.e.  0.1329   0.1071   0.0961   0.0746
## 
## sigma^2 estimated as 325309282:  log likelihood = -1890.8,  aic = 3789.6

checkresiduals(Propuesta.arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(0,1,1)[12]
## Q* = 14.447, df = 20, p-value = 0.8071
## 
## Model df: 4.   Total lags used: 24

Al parecer con la prueba Ljung-Box el valor de p-value es de 0.80, así que esta aceptando que no hay correlación en los residuos agrupados. Sin embargo la gráfica de los residuos muestra que aun hay un pico fuera de estacionalidad.

lo que comprueba que Auto.arima no siempre escoge el mejor modelo.

9.Presente la ecuación final y describa.

La mejor opción de modelado: Propuesta ARIMA(1,0,1)(1,1,1) #1

\[ Y_t=\phi Y_{t-1}- \theta Y_{t-1}+ \Phi Y_{t-12}- \Theta Y_{t-12}+ e_t \]

\[ Y_t=0.7428 Y_{t-1}- 0.2946 Y_{t-1}+ 0.0684 Y_{t-12}- 0.8273 Y_{t-12} \]

Concluyendo así que la mejor propuesta de modelado de la producción pesquera en México, es un ARIMA (1,0,1)(1,1,1). Mostrando que el modelo tiene una diferencia estacional y componentes estacionales y no estacionales AR(1) Y MA(1).