LA siguiente serie de tiempo contiene los resultados de los ??ndices de la productividad laboral nacional, proporcionados por la econom??a global (actividades primarias, secundarias y terciarias) toda en conjunto; estos datos en forma de indice representan datos existentes de producci??n, ventas, ingresos, empleo, horas trabajadas y remuneraciones provenientes de diversas fuentes estad??sticas del INEGI.
La forma medida de la base es en trimestres y la base de datos comienza en el primer trimestre del ano 2005 y finaliza en el primer trimestre del 2018. Podemos observar que tiene una tendencia positiva y una estacionalidad en el cuartos cuatrismestres de cada a??o donde por la temporada de dia de muertos y navidena la productividad laboral incremente. La serie sigue este `patron hasta el a??o 2008 y 2009, donde el desempleo estuvo fuerte por la crisis. Pero despues de este a??o para inicios de 2010 la serie se estabiliza y sigue creciente.
En esta grafica podemos observar que hay una tendencia positiva, exceptuando en el a??o 2009 donde hay una ca??da, pero la tendencia se sigue respetando a partir de ese a??o. Aparenta estacionalidad
library(ggplot2)
library(forecast)
library(seasonal)
library(fpp2)
## Loading required package: fma
## Loading required package: expsmooth
library(urca)
library(TSA)
## Loading required package: leaps
## Loading required package: locfit
## 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-23. For overview type 'help("mgcv-package")'.
## Loading required package: tseries
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
base<-read.csv('/Users/jesusraulmartinezpina/Desktop/Productividad.csv')
attach(base)
p.la<-ts(base$Pl,frequency = 4,start = c(2005,1))
autoplot(p.la,series = "Productividad Laboral")
ggseasonplot(p.la)
ggseasonplot(p.la, polar=TRUE)
fit<-decompose(p.la,type = "multiplicative")
autoplot(fit)
Por las graficas polares podemos observar que la estacionalidad se encuentra mayormente en el cuatrimestre numero 4 y podemos observar tambien patrones de estaciponalidad con una tendencia marcada.
BoxCox.ar(p.la)
Esta prueba no es muy utila para nostros. El optimo se encuentra en -1 y 0, no podemos aplicar logaritmos ya que la serie de tiempo esta medida en indice, pero buscaremos estabilizarla con diferencias estacionaria y no estacionaria.
dns.ds<-diff(diff(p.la),4)
autoplot(dns.ds)
Asi se comporta nuestra serie estacionalizada (aplicada una serie estacional y otra no estacional).
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.61
summary(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.61
Podemos ver una perdida de la tendencia y de una de la estacionalidad.
ggAcf(dns.ds)
Para la prueba ACF (Funcion de Auto correlacion Normal), se muestran varios rezagos que salen del limite permitido. El que m??s nos preocuopa es en el punto cuatro por lo que decimos que puede ser un Ma(4).
ggPacf(dns.ds)
Para la PACF (Prueba de Auto correlacion Parcial), se puede observar que hay 2 rezagos al principio que salen uno en el punto 4 y otro en el 7. Por lo que proponemos un Ar(4).
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 o
Para la EACF (Prueba de Auto correlacion Estacional) no est?? clara su interpretacon, los 0 estan muy dispersos. Probablemente probaremos con MA(4)y Ar(4) para nuestro primer modelo.
Observacion general
ggtsdisplay(dns.ds)
#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.3819 -0.4048 0.2388 0.1923 0.7433 -0.7926 -0.977 0.0286
## s.e. 1.8379 2.9337 1.1953 0.2715 1.8368 1.7500 1.568 1.9476
##
## sigma^2 estimated as 1.515: log likelihood=-83.77
## AIC=185.54 AICc=190.28 BIC=202.38
checkresiduals(fit1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0)(4,1,4)[4]
## Q* = 13.749, df = 3, p-value = 0.003268
##
## Model df: 8. Total lags used: 11
Para la propuesta 1 el modelo parece bueno , pero los rezagos en el puento 6 en la prueba ACF se sale del modelo y los residuos se salen. La prueba AIC = es de 185.54, esto es significativa, pero buscaremos mejor otro modelo.
# 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.84
checkresiduals(fit2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,0)[4]
## Q* = 8.5442, df = 5, p-value = 0.1287
##
## Model df: 3. Total lags used: 8
Para la propuesta 2 del modelo los rezagos son controlados en la prueba ACF, pero los residuos se descontrolan. La prueba AIC = es de 182.36, esto no es tan diferente al modelo anterior.
# 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.27
checkresiduals(fit3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(0,1,1)[4]
## Q* = 18.046, df = 6, p-value = 0.006118
##
## Model df: 2. Total lags used: 8
Podemos ver que los par??metros de la gr??fica ACF est??n controladas por excepci??n de peque??os puntos en la parte 5 y 8. Pero los residuos se controlaron mucho mejor. Aunque se puede ver que tal vez se deber??an controlar un poco m??s los rezagos. La prueba Aker aumento en 2 puntos, pero no es muy significativo en comparaci??n de las mejoras del modelo.
checkresiduals(fit3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(0,1,1)[4]
## Q* = 18.046, df = 6, p-value = 0.006118
##
## Model df: 2. Total lags used: 8
Ljung-Box test
data: Residuals from ARIMA(1,1,0)(0,1,1)[4] Q* = 18.046, df = 6, p-value = 0.006118
Model df: 2. Total lags used: 8
El valor nos da un valor p>0.05. Esto quiere decir que se no se rechaza la hipotesis nula y existe de correlacion serial. Se controlan mejor las puntos de la grafica ACP y hay un mejor control de los resiudos que en los casos anteriores
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.37
checkresiduals(fit4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(1,1,0)[4]
## Q* = 12.795, df = 6, p-value = 0.04641
##
## Model df: 2. Total lags used: 8
A comparaci??n del modelo Autorima en R studio, prefiero m??s nuestro modelo. La prueba p es mayor, pero los residuos son mas controlados en nuestro modelo. Por eso pensamos que es mejor nuestra opcion.
autoplot(rstandard(fit3))+ylab("Residuos Estandarizados")+geom_hline(yintercept = 0)
Como tal los residuos se comportan bien, pero necesitariamos ver los comportamientos de los a??os 2008 y 2009, donde hay una excepci??n y deberiamos buscar modelar especificamente esa excepcion.
qqnorm(residuals(fit3));qqline(residuals(fit3))
En el grafico quantil quantil, al principio podemos ver que se sale mas o menos de la media y al final de la serie que los puntos salen por debajo de la media.
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.
El pronostico de ARIMA, genera un pronostico a la media con tendencia mas regular.
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("Productividad laboral")+xlab("Year")+ylab("Productividad")+guides(colour=guide_legend(title="Forecast"))
library(seasonalview)
##
## Attaching package: 'seasonalview'
## The following object is masked from 'package:seasonal':
##
## view
base1<- ts(base, start = c(2005,1), frequency = 4)
plot(base1)
Vemos el avance de las variabes del indice del PIB y de la productividad laboral.
Se establecen comos series de tiempo.
p.laproductividad<- ts(base$Pl, start = c(2005,1), frequency = 4)
p.lapib<- ts(base$PIB, start = c(2005,1), frequency = 4)
Se puede ver que el indice del PIB tiene tendencia. Vamos a aplicar una diferencia
p.laproductividaddiff <- diff(diff(p.laproductividad, 4))
p.lapibdif <- diff(p.lapib)
Las juntamos en una base de datos
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
base2 <- cbind.zoo(Pl= p.laproductividaddiff, PIB = p.lapibdif)
plot(base2)
La serie de tiempo de la productividad laboral le falta el primer a??o, ya que se aplico las diferencias estacionales. Las ajustaremos reduciendo 4 trimestres.
base3 <- base2[-c(1:4),]
plot(base3)
Seguiremos con el pronostico de los rezagos
library(vars)
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
##
## cement, housing, petrol
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: lmtest
var1<-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.03472
El valor p es de 0.03472, no es muy significativo por lo que probaremos con 2 rezagos
var2<-VAR(base3,p=2,type="const")
serial.test(var2,lags.pt=10,type="PT.asymptotic")
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var2
## Chi-squared = 52.806, df = 32, p-value = 0.01174
El valor p es de 0.01174, no es muy significativo por lo que probaremos con 3 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.05162
Este rezago es el mejor, con 4 aumenta hasta en 0.4, pero eso quitaria muchos datos.
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.
Observamos los residuos del la Productividad laborla
plot(var3.serial,names="Pl")
cusum<-stability(var3,type="OLS-CUSUM")
plot(cusum)
rec.cusum<-stability(var3,type="Rec-CUSUM")
plot(rec.cusum)
causality(var3, cause='Pl')$Granger
##
## Granger causality H0: Pl do not Granger-cause PIB
##
## data: VAR object var3
## F-Test = 0.048243, df1 = 3, df2 = 76, p-value = 0.9859
causality(var3, cause='PIB')$Granger
##
## Granger causality H0: PIB do not Granger-cause Pl
##
## data: VAR object var3
## F-Test = 0.87591, df1 = 3, df2 = 76, p-value = 0.4575
No se reachaza la hipotesis nula de Granfer Causalidad. Por lo que el PIB esta relacionada con la Productividad laboral.
predictions<-predict(var3,n.ahead=10,ci=0.95)
class(predictions)
## [1] "varprd"
plot(predictions)
fanchart(predictions,names="Productividad laboral")
## Warning in fanchart(predictions, names = "Productividad laboral"):
## Invalid variable name(s) supplied, using first variable.
plot(irf(var3,impulse="Pl",response="PIB",ortho=T))
casi no se nota el impulso
plot(irf(var3,impulse="PIB",response="Pl",ortho=T))
aqui tiene un efecto mejor