TRABAJO FINAL DE SERIES DE TIEMPO

SERIES DE TIEMPO

Navarijo García Ricardo Daniel

Castillejos Arteaga Luis Felipe

Tasa de Desempleo en México

library(TSA) library(ggplot2) library(ggfortify) library(forecast) library(fpp2) library(forecast) library(fma) library(expsmooth) library(seasonal) library(urca) library(seasonal)

TABLA DE CONTENIDOS

1. Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc

2. Graficar la serie de tiempo, analizar patrones y observaciones atípicas. Apoye su análisis con una descomposición clásica.

3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza. Presentar grafica

4. Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadrática y/o medias estacionales. Presente el diagnóstico del modelo elegido.

5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales. Presentar grafica

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

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

8. Como parte del diagnóstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

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

10. Presente la ecuación final y describa.

11. Realizar el pronóstico ARIMA para dos años.

12. Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

1. Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc.

La tasa de desempleo, también conocida como tasa de paro, mide el nivel de desocupación en relación a la población activa. En otras palabras, es la parte de la población que estando en edad, condiciones y disposición de trabajar -población activa- no tiene puesto de trabajo.

La tasa de desempleo es muy útil para conocer las personas que no están trabajando. Su fórmula de cálculo es la población de 16 años y más que no está trabajando y busca trabajo, dividido entre la población económicamente activa de 16 años y más, esto es, ocupados más desocupados. Se calcula de la siguiente manera: tasa de desempleo=(No. de desempleados/Poblacion activa)*100 datos sacdos de http://www.fao.org

2. Graficar la serie de tiempo, analizar patrones y observaciones atípicas. Apoye su análisis con una descomposición clásica.

base<-read.csv(file.choose()) tdd<-ts(base$tdd, start= c(1987,1),frequency= 12)

library(ggplot2)
library(ggfortify)
autoplot(tdd)

aquí podemos apreciar como tenemos una deserie donde aument eldesempleo en los 90’s por el efecto tequila, donde tuvimos una depresicación de la moneda mexicana, de igual manera podemos observar como baja poco a poco hasta el año 2000 donde tenemos un incremento hasta llegar a otra crisis en el años 2008-2009 pero no tan fuerte como en los 90’s

fit <- decompose(tdd, type=‘multiplicative’)

Para el caso clásico obtenemos los componentes separados de nuestra serie, en relacion a la estacionalidad, tendencia y residuos.

fit2 <- seas(tdd, x11=“”)

Se realizo una gráfica que compara el comportamiento de nuestros datos con el de un ajuste estacional y observamos un comportamiento más preciso en el componente tendencial, además de que la estacionalidad refleja cambios que van acortando los extremos verticales con una pequeña alza a mitad del periodo, que después regresa a su actividad inicial.

3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza. Presentar grafica

En nuestros caso se estan trabajando con índices, por lo cual no es correcto el hacer uso de logaritmos, por lo que solo consideramos la aplicación de la diferenciación. Por lo anterior dicho, no hacemos consideración de la transformación Box Cox.

4. Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadrática y/o medias estacionales. Presente el diagnóstico del modelo elegido.

modelo1 <- lm(tdd~time(tdd))

summary(modelo1)

Call:
lm(formula = tdd ~ time(tdd))

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6444 -0.6376 -0.1710  0.5324  4.0430 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -79.562984  10.827546  -7.348 1.28e-12 ***
time(tdd)     0.041652   0.005407   7.704 1.20e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9445 on 373 degrees of freedom
Multiple R-squared:  0.1373,    Adjusted R-squared:  0.135 
F-statistic: 59.35 on 1 and 373 DF,  p-value: 1.205e-13

autoplot(modelo1, which = 1:6, ncol = 2, label.size = 3)

package <U+393C><U+3E31>bindrcpp<U+393C><U+3E32> was built under R version 3.4.4

shapiro.test(rstudent(modelo1))

    Shapiro-Wilk normality test

data:  rstudent(modelo1)
W = 0.93097, p-value = 3.725e-12

mes. <- season(tdd) modelo2 <- lm(tdd~mes.)

summary(modelo2)

Call:
lm(formula = tdd ~ mes.)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7621 -0.7700 -0.2213  0.7690  3.6787 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    3.827749   0.182065  21.024   <2e-16 ***
mes.February   0.034362   0.257479   0.133    0.894    
mes.March     -0.005961   0.257479  -0.023    0.982    
mes.April      0.050283   0.259547   0.194    0.846    
mes.May        0.050119   0.259547   0.193    0.847    
mes.June       0.046411   0.259547   0.179    0.858    
mes.July       0.062057   0.259547   0.239    0.811    
mes.August     0.093549   0.259547   0.360    0.719    
mes.September  0.002964   0.259547   0.011    0.991    
mes.October   -0.015663   0.259547  -0.060    0.952    
mes.November   0.003626   0.259547   0.014    0.989    
mes.December  -0.071200   0.259547  -0.274    0.784    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.03 on 363 degrees of freedom
Multiple R-squared:  0.001701,  Adjusted R-squared:  -0.02855 
F-statistic: 0.05622 on 11 and 363 DF,  p-value: 1

shapiro.test(rstudent(modelo2))

    Shapiro-Wilk normality test

data:  rstudent(modelo2)
W = 0.95154, p-value = 9.308e-10

modelo3 <- lm(tdd~mes. -1)

summary(modelo3)

Call:
lm(formula = tdd ~ mes. - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7621 -0.7700 -0.2213  0.7690  3.6787 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
mes.January     3.8277     0.1821   21.02   <2e-16 ***
mes.February    3.8621     0.1821   21.21   <2e-16 ***
mes.March       3.8218     0.1821   20.99   <2e-16 ***
mes.April       3.8780     0.1850   20.96   <2e-16 ***
mes.May         3.8779     0.1850   20.96   <2e-16 ***
mes.June        3.8742     0.1850   20.94   <2e-16 ***
mes.July        3.8898     0.1850   21.03   <2e-16 ***
mes.August      3.9213     0.1850   21.20   <2e-16 ***
mes.September   3.8307     0.1850   20.71   <2e-16 ***
mes.October     3.8121     0.1850   20.61   <2e-16 ***
mes.November    3.8314     0.1850   20.71   <2e-16 ***
mes.December    3.7565     0.1850   20.31   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.03 on 363 degrees of freedom
Multiple R-squared:  0.9352,    Adjusted R-squared:  0.933 
F-statistic: 436.4 on 12 and 363 DF,  p-value: < 2.2e-16

shapiro.test(rstudent(modelo3))

    Shapiro-Wilk normality test

data:  rstudent(modelo3)
W = 0.95154, p-value = 9.308e-10

5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales. Presentar grafica

tdd.d<-autoplot(diff(tdd))

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

Diferenciar hasta que la serie sea estacionaria. Ho: Existe raíz unitaria Ha: No existe raíz unitaria

adf.test(tdd)

    Augmented Dickey-Fuller Test

data:  tdd
Dickey-Fuller = -2.5049, Lag order = 7, p-value = 0.3639
alternative hypothesis: stationary

adf.test(tdd.d)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  tdd.d
Dickey-Fuller = -6.0305, Lag order = 7, p-value = 0.01
alternative hypothesis: stationary

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

eacf(tdd.d)

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

8. Como parte del diagnóstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

fit1<-Arima(tdd.d, order=c(1,1,0), seasonal = c(1,1,0))

fit1

Series: tdd.d 
ARIMA(1,1,0)(1,1,0)[12] 

Coefficients:
          ar1     sar1
      -0.6851  -0.5895
s.e.   0.0388   0.0419

sigma^2 estimated as 0.1467:  log likelihood=-167.61
AIC=341.23   AICc=341.29   BIC=352.89

checkresiduals(fit1)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(1,1,0)[12]
Q* = 197.66, df = 22, p-value < 2.2e-16

Model df: 2.   Total lags used: 24

fit2<-Arima(tdd.d, order=c(2,1,0), seasonal = c(1,1,0))

fit2

Series: tdd.d 
ARIMA(2,1,0)(1,1,0)[12] 

Coefficients:
          ar1      ar2     sar1
      -0.9509  -0.3918  -0.5877
s.e.   0.0488   0.0485   0.0423

sigma^2 estimated as 0.1245:  log likelihood=-137.73
AIC=283.46   AICc=283.57   BIC=299.01

checkresiduals(fit2)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,0)(1,1,0)[12]
Q* = 153.02, df = 21, p-value < 2.2e-16

Model df: 3.   Total lags used: 24

fit3<-auto.arima(tdd.d)

fit3

Series: tdd.d 
ARIMA(3,0,5)(0,0,2)[12] with zero mean 

Coefficients:
          ar1     ar2     ar3     ma1      ma2      ma3     ma4     ma5     sma1    sma2
      -1.4183  0.1019  0.5373  1.1735  -0.3911  -0.3383  0.3859  0.1442  -0.1353  0.1102
s.e.   0.2505  0.4935  0.2457  0.2563   0.4246   0.1610  0.0761  0.0689   0.0546  0.0555

sigma^2 estimated as 0.05748:  log likelihood=7.77
AIC=6.46   AICc=7.19   BIC=49.63

checkresiduals(fit3)

    Ljung-Box test

data:  Residuals from ARIMA(3,0,5)(0,0,2)[12] with zero mean
Q* = 13.83, df = 14, p-value = 0.4625

Model df: 10.   Total lags used: 24

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

fit4<-auto.arima(tdd.d, stepwise = FALSE, approximation = FALSE)

fit4

Series: tdd.d 
ARIMA(1,0,2)(1,0,0)[12] with zero mean 

Coefficients:
         ar1      ma1     ma2     sar1
      0.7332  -1.0379  0.3424  -0.1377
s.e.  0.1178   0.1166  0.0512   0.0521

sigma^2 estimated as 0.05974:  log likelihood=-1.96
AIC=13.93   AICc=14.09   BIC=33.55

checkresiduals(fit4)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,2)(1,0,0)[12] with zero mean
Q* = 23.12, df = 20, p-value = 0.2829

Model df: 4.   Total lags used: 24

10. Presente la ecuación final y describa.

fit5<-Arima(tdd.d, order=c(3,0,5), seasonal = c(0,0,2))

fit5

Series: tdd.d 
ARIMA(3,0,5)(0,0,2)[12] with non-zero mean 

Coefficients:
          ar1     ar2    ar3     ma1      ma2      ma3     ma4     ma5     sma1    sma2
      -1.4143  0.1106  0.542  1.1685  -0.3997  -0.3421  0.3842  0.1426  -0.1356  0.1089
s.e.   0.2425  0.4779  0.238  0.2483   0.4105   0.1565  0.0762  0.0686   0.0545  0.0555
         mean
      -0.0050
s.e.   0.0131

sigma^2 estimated as 0.05762:  log likelihood=7.84
AIC=8.32   AICc=9.18   BIC=55.41

checkresiduals(fit5)

    Ljung-Box test

data:  Residuals from ARIMA(3,0,5)(0,0,2)[12] with non-zero mean
Q* = 13.806, df = 13, p-value = 0.3876

Model df: 11.   Total lags used: 24

Se ecoge este modelo porque es donde nuestro ACF no se salen y los residuos no se salen mucho a comparación de otros

11. Realizar el pronóstico ARIMA para dos años.

getrmse<-function(x,h,…) { train.end<-time(x)[length(x)-h] test.start<-time(x)[length(x)-h+1] train<-window(x,end=train.end) test<-window(x,start=test.start) fit<-Arima(train,…) fc<-forecast(fit,h=h) return(accuracy(fc,test)[2,“RMSE”]) }

getrmse(tdd,h=24,order=c(1,0,2),seasonal=c(1,0,0))

[1] 0.4559535

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,0))

[1] 0.617638

getrmse(tdd,h=24,order=c(2,1,0),seasonal=c(1,1,0))

[1] 0.4229213

getrmse(tdd,h=24,order=c(2,1,0),seasonal=c(1,1,2))

[1] 0.4072446

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,2))

[1] 0.2503558

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,1))

[1] 0.2523106

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(0,1,1))

[1] 0.2550056

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(1,1,0))

[1] 0.3094996

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(1,1,1))

[1] 0.3900453

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(0,1,2))

[1] 0.4088291

getrmse(tdd,h=24,order=c(3,1,4),seasonal=c(1,0,0))

[1] 0.3602454

getrmse(tdd,h=24,order=c(3,0,5),seasonal=c(0,0,2))

[1] 0.2359293

getrmse(tdd,h=24,order=c(3,1,1),seasonal=c(0,0,2))

[1] 0.3267824

fit6<-Arima(tdd, order=c(2,1,1),seasonal=c(1,1,2))

fit6

Series: tdd 
ARIMA(2,1,1)(1,1,2)[12] 

Coefficients:
         ar1     ar2      ma1     sar1     sma1     sma2
      0.5434  0.3153  -0.8153  -0.2883  -0.8138  -0.1276
s.e.  0.0999  0.0517   0.0927   0.3488   0.3562   0.3553

sigma^2 estimated as 0.06294:  log likelihood=-25.08
AIC=64.16   AICc=64.48   BIC=91.4

autoplot(forecast(fit6,h=12,type='b',xlab='Time',ylab='Color  Property'))

The non-existent type arguments will be ignored.The non-existent xlab arguments will be ignored.The non-existent ylab arguments will be ignored.

12. Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

Para poder realizar el pronostico de nuestra serie, primeramente se realizo un recorte de los ultimos dos años para despues pronosticar nuestra serie mediante los metodos mean(media), naive y naive estacional.

tdd.c <- window(tdd,start=c(1987,4),end=c(2015,12))

http://rpubs.com/Richie1995/ http://rpubs.com/Richie1995/397935