Trabajo Final(Productividad laboral)
Martinez Piña Jesus Raul
Introducción
En ésta ocasión presentar una serie de tiempo en donde se analiza mediante tasas de desocupación (mostradas por INEGI) el nivel de desempleo que ha habido en M?xico para las personas con preparatoria y universidad terminada.
El periodo es anual y la base de datos comienza en enero del a?o 2005 y finaliza en marzo del 2018 Como an?lisis Econ?mico hay que recordar que durante ese año, fue el último en donde nuestro Gobierno estaba bajo el sexenio de Vicente Fox, y al llegar Felipe Calderón una de sus propuesta es que el desempleo iba a descender. Por lo que se espera que nuestra gráfica tenga una tendencia negativa. Al igual veremos cómo fue su comportamiento durante el sexenio de Enrique Peña Nieto
Comportamiento de la Serie
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4library(forecast)## Warning: package 'forecast' was built under R version 3.4.4library(seasonal)## Warning: package 'seasonal' was built under R version 3.4.4library(fpp2)## Warning: package 'fpp2' was built under R version 3.4.4## Loading required package: fma## Warning: package 'fma' was built under R version 3.4.4## Loading required package: expsmooth## Warning: package 'expsmooth' was built under R version 3.4.4library(urca)## Warning: package 'urca' was built under R version 3.4.4library(TSA)## Warning: package 'TSA' was built under R version 3.4.4## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.4.4## Loading required package: locfit## Warning: package 'locfit' was built under R version 3.4.4## locfit 1.5-9.1 2013-03-22## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:forecast':
##
## getResponse## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.## Loading required package: tseries## Warning: package 'tseries' was built under R version 3.4.4##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarbase<-read.csv('C:/Users/UDI/Documents/Productividad.csv')
attach(base)
p.la<-ts(base$Pl,frequency = 4,start = c(2005,1))
autoplot(p.la,series = "Productividad Laboral")
Claramente se observa una tendencia positiva de nuestra serie que se ve más marcada a partir del 2010 en adelante. Ahora haremos una descomposici?n a nuestra Serie.
Análisis y Descomposición
fit<-decompose(p.la,type = "multiplicative")
autoplot(fit)
Hay una clara tendencia en la base de datos hacia el final de nuestra serie, al principio no es muy marcado, pero llega un punto en donde s?lo es ascendente. Adem?s de mostrar ciertos patrones de estacionalidad.
Prueba BoxCox
BoxCox.ar(p.la)
El ?ptimo se encuentra en 1, as? que aplicaremos logaritmos para estabilizar las varianzas
No aplicamos logaritmos, ya que la serie de timpo es en base a un indice.
Prueba Dickey Fuller
summary(ur.df(p.la,type = "trend",selectlags = "AIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8935 -1.2068 -0.2877 1.2775 4.1386
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.27628 14.45007 2.926 0.00528 **
## z.lag.1 -0.43528 0.14907 -2.920 0.00536 **
## tt 0.04710 0.02457 1.917 0.06134 .
## z.diff.lag -0.20191 0.14143 -1.428 0.16000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.958 on 47 degrees of freedom
## Multiple R-squared: 0.3051, Adjusted R-squared: 0.2607
## F-statistic: 6.877 on 3 and 47 DF, p-value: 0.000619
##
##
## Value of test-statistic is: -2.92 2.9104 4.2634
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47En ?sta, la serie ?nicamente tiene logaritmos, pero sigue presentando tendencias. Nos presenta 3 valores cr?ticos y phis. Por lo que decid? quitar la tendencia de la serie usando diferenciaci?n.
Aplicando Diferencias Estacionales y con Tendencia
dns.ds<-diff(diff(p.la),4)
autoplot(dns.ds)
As? luce nuestra serie una vez que aplicamos diferencias y logaritmos, como podemos ver, ya qued? estacionalizada.
summary(ur.df(dns.ds,type="trend",selectlags = "AIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8510 -0.9920 -0.1483 1.2686 6.1815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.129306 0.629308 -0.205 0.838
## z.lag.1 -1.068259 0.234622 -4.553 4.48e-05 ***
## tt 0.002941 0.022608 0.130 0.897
## z.diff.lag -0.053970 0.165404 -0.326 0.746
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.03 on 42 degrees of freedom
## Multiple R-squared: 0.5501, Adjusted R-squared: 0.518
## F-statistic: 17.12 on 3 and 42 DF, p-value: 2.061e-07
##
##
## Value of test-statistic is: -4.5531 6.9178 10.3761
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.15 -3.50 -3.18
## phi2 7.02 5.13 4.31
## phi3 9.31 6.73 5.61summary(ur.df(dns.ds,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
## -4.947 -1.061 -0.215 1.180 6.093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -1.06613 0.22909 -4.654 3e-05 ***
## z.diff.lag -0.05568 0.16148 -0.345 0.732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.985 on 44 degrees of freedom
## Multiple R-squared: 0.5497, Adjusted R-squared: 0.5292
## F-statistic: 26.85 on 2 and 44 DF, p-value: 2.384e-08
##
##
## Value of test-statistic is: -4.6538
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.62 -1.95 -1.61Al momento de hacerlo, identificamos que la serie pierde su tendencia (como se mostr? en la gr?fica anterior). Adem?s, que ?nicamente presenta 1 valor cr?tico para 3 valores de Tau, que era lo que se buscaba inicialmente. Por lo que la serie diferenciada es la elegida para la aplicar las pruebas de correlaci?n.
Prueba ACF
ggAcf(dns.ds)
Para la prueba ACF (Funci?n de Auto correlaci?n Normal), claramente muestra varios rezagos que salen del l?mite permitido. Por lo que el modelo elegido para ?ste caso es un MA2.
Prueba PACF
ggPacf(dns.ds)
Para la PACF (Prueba de Auto correlaci?n Parcial), claramente se ve que hay 3 rezagos al principio y en medio de las muestras que salen del l?mite permitido. Por lo que el modelo elegido para ?ste caso es una AR3
Prueba EACF
eacf(dns.ds)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o x o o o x o o o o o o
## 1 x o o x o o o o o o o o o o
## 2 o o o x o o o o o o o o o o
## 3 o o o x o o o o o o o o o o
## 4 o o x o o o o o o o o o o o
## 5 x o x o o o o o o o o o o o
## 6 x o o x o o o o o o o o o o
## 7 o o o o o o o o o o o o o oPara la EACF (Prueba de Auto correlaci?n Estacional) es algo confusa su interpretaci?n, ya que los 0 est?n muy dispersos. Por lo tanto, probaremos con los modelos MA2 y AR2
ggtsdisplay(dns.ds)
Aqu? una mejor observaci?n de todas las gr?ficas juntas
Propuesta 1
#Primer modelo propuesto ARIMA (0,1,0) (4,1,4)
fit1 <- Arima( p.la,order=c(0,1,0),seasonal=c(4,1,4))
fit1## Series: p.la
## ARIMA(0,1,0)(4,1,4)[4]
##
## Coefficients:
## sar1 sar2 sar3 sar4 sma1 sma2 sma3 sma4
## -1.3853 -0.4119 0.2351 0.192 0.7449 -0.7889 -0.9789 0.0248
## s.e. 1.8452 2.9451 1.2028 0.271 1.8490 1.7555 1.5762 1.9538
##
## sigma^2 estimated as 1.516: log likelihood=-83.77
## AIC=185.54 AICc=190.28 BIC=202.38checkresiduals(fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(4,1,4)[4]
## Q* = 13.767, df = 3, p-value = 0.00324
##
## Model df: 8. Total lags used: 11Para el modelo 1, al cual llamaremos “Propuesta 1” iniciamos aplicando la funci?n “arima” creando as? el modelo AR2, el cual nos arroja los siguientes datos: . Si dividimos -.51/.07, nos da una relaci?n bastante baja . El resultado de AIC es de -333.26
Propuesta 2
# Segundo modelo propuesto ARIMA (1,1,1)(1,1,0)
fit2 <- Arima( p.la,order=c(1,1,1),seasonal=c(1,1,0))
fit2## Series: p.la
## ARIMA(1,1,1)(1,1,0)[4]
##
## Coefficients:
## ar1 ma1 sar1
## -0.6950 1.0000 -0.6283
## s.e. 0.1312 0.0616 0.1157
##
## sigma^2 estimated as 2.132: log likelihood=-87.18
## AIC=182.36 AICc=183.29 BIC=189.84Ahora para el siguiente modelo, al cual llamaremos propuesta 2, iniciamos aplicando la funci?n “arima” indicando ahora que queremos un MA2. El cual no da los siguientes datos: . Una correlaci?n en MA2 mayor a la que ten?amos anteriormente . Un AIC mayor al anterior
Por lo que podemos afirmar que el mejor modelo para ?sta serie de tiempo es un: AR2
# Tercer modelo propuesto
fit3 <- Arima( p.la,order=c(1,1,0),seasonal=c(0,1,1))
fit3## Series: p.la
## ARIMA(1,1,0)(0,1,1)[4]
##
## Coefficients:
## ar1 sma1
## -0.0817 -0.9998
## s.e. 0.1511 0.5905
##
## sigma^2 estimated as 1.999: log likelihood=-88.83
## AIC=183.66 AICc=184.2 BIC=189.27Prubea Ljung-Box ( del mejor)
checkresiduals(fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(4,1,4)[4]
## Q* = 13.767, df = 3, p-value = 0.00324
##
## Model df: 8. Total lags used: 11Ljung-Box test
data: Residuals from ARIMA(2,0,0) with non-zero mean Q* = 91.005, df = 21, p-value = 1.082e-10
Model df: 3. Total lags used: 24
Al revisar los residuos de nuestro modelo AR2, nos da un valor p<0.05. Esto quiere decir que se acepta la hip?tesis nula y por lo tanto la existencia de correlaci?n serial
Si bien no es perfecta, muestra de una manera bastante aceptable en donde el histograma se ve de buena forma a trav?s de la campana. Los residuos son consistententes en sus varianzas, aunque a?n existe un rezago que sale de los l?mites permitidos.
Propuesta Auto() Arima
fit4<-auto.arima(p.la)
fit4## Series: p.la
## ARIMA(1,0,0)(1,1,0)[4]
##
## Coefficients:
## ar1 sar1
## 0.8127 -0.5934
## s.e. 0.0843 0.1197
##
## sigma^2 estimated as 2.187: log likelihood=-88.85
## AIC=183.7 AICc=184.23 BIC=189.37Para Rstudio, el mejor modelo a elegir es un ARMA (2,0,1) As? que hicimos una buena elecci?n en el modelo. Ahora veremos el comportamiento de los residuos en diferentes gr?ficas.
Comportamiento de los Residuos
autoplot(rstandard(fit3))+ylab("Residuos Estandarizados")+geom_hline(yintercept = 0)
Al principio de la serie hay algunos puntos que exceden inusualmente del n?mero 2. Habr?a que ver si hay otros factores que afecten el desempleo de las personas con Universidad terminada, con condiciones sociales durante 2005-2008
Gr?fico Quantil-Quantil
qqnorm(residuals(fit3));qqline(residuals(fit3))
En el gr?fico quantil quantil, al principio y al final de la serie se ve claramente los puntos que salen (por muy poco) de la media.
Pron?sticos ARIMA
autoplot(forecast(fit4,h=4,type="b",xlab="Time",ylab="Desempleo"))## Warning in forecast.Arima(fit4, h = 4, type = "b", xlab = "Time", ylab
## = "Desempleo"): The non-existent type arguments will be ignored.The non-
## existent xlab arguments will be ignored.The non-existent ylab arguments
## will be ignored.
Como podemos ver, en el pron?stico de ARIMA, lo ?nico que hace es generar un pronostico que tienda a la media hasta hacer su movimiento m?s regular.
Prono?stico Simple
autoplot(p.la)+forecast::autolayer(meanf(p.la,h=3),PI=FALSE,series="Mean")+forecast::autolayer(naive(p.la,h=3),PI=FALSE,series="Na?ve")+forecast::autolayer(snaive(p.la,h=3),PI=FALSE,series="Seasonal na?ve")+ggtitle("Desempleo Juvenil")+xlab("Year")+ylab("Desempleo")+guides(colour=guide_legend(title="Forecast"))
Teor?a Keynesiana
El objetivo de poner la Teor?a Keynesiana es porque dentro de ella se contempla un factor que se llama Gasto de Gobierno, El gasto p?blico contempla todas las erogaciones que realiza el Estado. El gasto p?blico puede clasificarse en funci?n a la naturaleza de los servicios que las instituciones p?blicas brindan a la comunidad. Los gastos clasificados por finalidad permiten determinar los objetivos generales y las acciones a trav?s de las cuales se estima alcanzarlos.
Dicho de otro modo; para Keynes, el Gobierno es el agente que regula la econom?a, incetiva su partici?ci?n dentro de la misma, adem?s mandejaba un concepto que denomin? el pleno empleo. Al relacionarlos, Keynes dec?a que entre m?s participara el gobierno, menos desempleo existir?a, pero entre menos gastara, m?s desempleo habr?a.
Nueva Variable Inclu?da
library(seasonalview)## Warning: package 'seasonalview' was built under R version 3.4.4##
## Attaching package: 'seasonalview'## The following object is masked from 'package:seasonal':
##
## viewp.la1<- ts(base, start = c(2005,1), frequency = 4)
plot(p.la1)
Aqui podemos ver el avance de ambas variables (Desempleo Juvenil y Gasto de Gobierno).
Ahora declaramos cada variable como serie de tiempo.
p.laproductividad<- ts(base$Pl, start = c(2005,1), frequency = 4)
p.lapib<- ts(base$PIB, start = c(2005,1), frequency = 4)Estabilizaci?n de las Variables
Ahora debemos diferenciar para cada variable lo que se necesite para estabilizarla
p.laproductividaddiff <- diff(diff(p.laproductividad, 4))
p.lapibdif <- diff(p.lapib)Ahora debemos guardar ambas variables en una sola base de Datos
library(zoo)## Warning: package 'zoo' was built under R version 3.4.4##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbase2 <- cbind.zoo(D= p.laproductividaddiff, G = p.lapibdif)
plot(base2)
Como podemos ver en la gr?fica, nuestra serie de tiempo de Desempleo est? cortada al principio, ya que como aplicamos diferencias estacionales, nos hacen falta 12 datos, por lo que debemos ajustarlas
base3 <- base2[-c(1:4),]
plot(base3)
Perfecto ambas variables est?n en la misma longitud por lo que podemos comenzar a pronosticar los rezagos.
Rezagos
Con 1 Rezago
library(vars)## Warning: package 'vars' was built under R version 3.4.4## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:fma':
##
## cement, housing, petrol## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 3.4.4## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 3.4.4## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 3.4.4var1<-VAR(base3,p=1,type="const")
serial.test(var1,lags.pt=10,type="PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var1
## Chi-squared = 52.843, df = 36, p-value = 0.03472Nuestro valor p es de 0.00013, es muy peque?o, por lo que probaremos con 2 rezagos
Con 2 Rezagos
var3<-VAR(base3,p=3,type="const")
serial.test(var3,lags.pt=10,type="PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var3
## Chi-squared = 41.187, df = 28, p-value = 0.05162Nuestro valor p es 0.20. Mejor? muchisimo a comparaci?n del anterior. As? que nos quedamos con 2 rezagos.
Prueba de Diagn?stico
var3.serial<-serial.test(var3,lags.pt=10,type="PT.asymptotic")
plot(var3.serial,names="Consumption")## Warning in plot.varcheck(var3.serial, names = "Consumption"):
## Invalid residual name(s) supplied, using residuals of first variable.
Ahora veremos los residuos del Desempleo
plot(var3.serial,names="Desempleo")## Warning in plot.varcheck(var3.serial, names = "Desempleo"):
## Invalid residual name(s) supplied, using residuals of first variable.
cusum<-stability(var3,type="OLS-CUSUM")
plot(cusum)
rec.cusum<-stability(var3,type="Rec-CUSUM")
plot(rec.cusum)
Prueba de Causalidad de Granger
causality(var3, cause='D')$Granger##
## Granger causality H0: D do not Granger-cause G
##
## data: VAR object var3
## F-Test = 0.048243, df1 = 3, df2 = 76, p-value = 0.9859causality(var3, cause='G')$Granger##
## Granger causality H0: G do not Granger-cause D
##
## data: VAR object var3
## F-Test = 0.87591, df1 = 3, df2 = 76, p-value = 0.4575Claramente no es posible rechazar la hip?tesis nula de no Granfer Causalidad. Por lo que el Gasto de Gobierno S? causa al Desempleo
Predicciones VAR
predictions<-predict(var3,n.ahead=10,ci=0.95)
class(predictions)## [1] "varprd"plot(predictions)
fanchart(predictions,names="Desempleo")## Warning in fanchart(predictions, names = "Desempleo"):
## Invalid variable name(s) supplied, using first variable.
Impulso - Respuesta
plot(irf(var3,impulse="D",response="G",ortho=T))
plot(irf(var3,impulse="G",response="D",ortho=T))
Como podemos ver en la primer gr?fica, alcanza valores de hasta 2.