Evidencia_a01284354

Econometría 2

Grupo 102

EVIDENCIA PRONÓSTICOS DE LA ECONOMÍA

INTRODUCCIÓN

De acuerdo con la Organización Internacional del Trabajo (OIT), la subocupación es definida como la situación en la cual los trabajadores reducen sus horas de trabajo y laboran en empleos menos satisfactorios con el objetivo de evitar el desempleo (Morales, s.f.). De acuerdo con un estudio realizado por Clara Judisman, quién fungió como Directora General del Registro Nacional de Electores en México, en las naciones de menor desarrollo resulta más efectivo el análisis de la subocupación en lugar de el análisis del desempleo, esto sucede debido a que de acuerdo con B. Datta el problema principal en estos países “no es el volumen, número o proporción de personas empleadas, sino el grado del empleo” (Judisman,s.f, p.269). En México, gran porcentaje de la población carece de recursos financieros para poder laborar en un trabajo adecuado, o con un periodo largo, por lo tanto, los individuos prefieren tomar cualquier empleo antes de permanecer desempleados (Judisman,s.f,p. 275). Por otro lado, Bantegui plantea que “la inadecuada utilización de la fuerza de trabajo no se expresa por el desempleo, sino, en la carencia de empleos satisfactorios” (Judisman,s.f,p. 275). Debido a esto, en el presente trabajo se realizará un pronóstico de la variable de subocupación con el programa R Studio. En este caso se utilizará una base de datos obtenida del INEGI, la cual cuenta con 70 observaciones. Por otro lado, la serie es de la variación trimestral de la tasa de subocupación en México, y se tomaron datos desde el primer trimestre del año 2005, al segundo trimestre del año 2022.

DESARROLLO

1. Modelar la serie de Tiempo

• Representación Gráfica
• Prueba Dickey Fuller
• Identificar el Modelo
• Validación de Modelo
• Refinación del Modelo

2. Pronóstico de la Tasa de Subocupación

• Estimación Auto Arima
• Generación de Grupos de Entrenamiento y Prueba
• Comparación de Modelos - Prueba Diabold Mariano
• Pronóstico

1. Modelar la Serie de Tiempo

Lectura de base de datos

datos <- read_excel("~/Desktop/MÉTODOS_ R/subocupacion.xlsx")
date = as.Date(datos$Fecha, format = "%Y/%m/%d")
subocupacion = datos$sub

Representación Gráfica

En primera instancia se realiza una representación gráfica para observar si la serie es estacionaria.

plot(datos$sub, type = "l", xlab = "", ylab = "%", main = "")

Sin embargo, el análisis gráfico no es suficiente para observar esta propiedad. Por lo tanto, se procede a realizar una prueba que contribuya a demostrar estacionariedad.

Prueba Dicky Fuller (nivel de significancia 5%)

El objetivo de esta prueba es observar si la serie a utilizar es estacionaria o no, esto quiere decir que los momentos estadísticos no dependan del tiempo. Esta prueba aumentada de Dickey-Fuller se estima con una regresión auxiliar. La ecuación asociada es:

\[\begin{align*} \Delta y_t = \beta_0 + \delta t + \phi y_{t-1} + \sum_{i=1}^{p}\beta \Delta y_{t-i}+e_t \end{align*}\]

La hipotesis asociada a la prueba es la siguiente: \[\begin{align*} H_0: \phi = 0 \\ H_a: \phi \neq 0 \end{align*}\]

Si se rechaza H0 la serie no tiene raíz unitaria, por lo tanto, la serie es estacionaria.

Se usará un nivel de significancia del 5%

Realizando por primera vez la prueba se observa lo siguiente:

sub.ts <- ts(datos ,start = c(2005,01), end = c(2022,02), frequency = 4)#Convirtiendo a series de tiempo
adf.sub.trend <- ur.df(sub.ts[,"sub"], type = "trend", selectlags = "BIC")
summary(adf.sub.trend)

############################################### 
# 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 
-1.6583 -0.7816 -0.2154  0.1732 15.8679 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.48845    0.91843   2.709  0.00864 **
z.lag.1     -0.37481    0.11449  -3.274  0.00171 **
tt           0.02132    0.01487   1.433  0.15670   
z.diff.lag  -0.11018    0.12419  -0.887  0.37829   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.143 on 64 degrees of freedom
Multiple R-squared:  0.219, Adjusted R-squared:  0.1824 
F-statistic: 5.981 on 3 and 64 DF,  p-value: 0.001168


Value of test-statistic is: -3.2736 3.5754 5.3605 

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.47

Debido a que el valor de test-statistic no es más negativo que el valor de tau3 en un nivel de significancia al 5%, NO se rechaza H0. Por otro lado, ya que el valor test-statistic 3.57 no es mayor al valor de phi2, en un nivel de significancia de 5%, NO se rechaza H0.

Debido a que la serie sigo contando con raíz unitaria, se realiza otra prueba removiendo la tendencia con la función “drift”.

adf.sub.drift <- ur.df(sub.ts[,"sub"], type = "drift", selectlags = "BIC")
summary(adf.sub.drift)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3886 -0.5929 -0.2755 -0.0318 16.4390 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   2.6016     0.9224   2.820  0.00635 **
z.lag.1      -0.3001     0.1027  -2.920  0.00480 **
z.diff.lag   -0.1464     0.1226  -1.195  0.23657   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.16 on 65 degrees of freedom
Multiple R-squared:  0.1939,    Adjusted R-squared:  0.1691 
F-statistic: 7.818 on 2 and 65 DF,  p-value: 0.000907


Value of test-statistic is: -2.9204 4.267 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1  6.70  4.71  3.86

Ya que el valor de test-statistic (-2.89) es más negativo que el valor de tau (-2.89) en un nivel de significancia al 5%, se rechaza H0.

Se concluye que la serie es estacionaria y NO tiene raíz unitaria. El modelo es estacionario y tiene un orden de integración de cero.

Identificar el Modelo

Correlograma

Con el objetivo de identificar el modelo adecuadamente se utilizará un correlograma, el cual permitirá observar si se cuenta con un modelo AR, MA o ARMA.

par(mfrow = c(1,2))
acf(subocupacion, lag = 20)
pacf(subocupacion, lag = 20)

Tomando en cuenta que el decaimiento geométrico se encuentra en ACF y ACP, se concluye que es un modelo ARMA (1,1). Por lo tanto, se estimará la función arima como un modelo de orden (1,0,1).

Validación Modelo ARMA (1,1)

Estimación Modelo ARMA(1,1)

modelo.arma <- arima(datos$sub, order = c(1, 0, 1)) 
modelo.arma

Call:
arima(x = datos$sub, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.7543  -0.1997     8.6170
s.e.  0.1099   0.1619     0.7814

sigma^2 estimated as 4.354:  log likelihood = -151.1,  aic = 310.19

summary(modelo.arma)

Call:
arima(x = datos$sub, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.7543  -0.1997     8.6170
s.e.  0.1099   0.1619     0.7814

sigma^2 estimated as 4.354:  log likelihood = -151.1,  aic = 310.19

Training set error measures:
                      ME     RMSE       MAE       MPE     MAPE      MASE
Training set 0.003161995 2.086719 0.7425833 -2.635958 7.224152 0.8729572
                     ACF1
Training set 0.0003464964

Prueba Individual

par(mfrow = c(1,2))
acf(modelo.arma$residuals, lag = 20)
pacf(modelo.arma$residuals, lag = 20)

Debido a que la estructura en el AC y ACP se encuentran dentro de las bandas de confianza, se concluye que es un modelo adecuado.

Prueba Conjunta

Por otro lado, en la Prueba Conjunta de Lyung, se busca observar si hay estructura de autocorrelación.

Box.test(modelo.arma$residuals, type = "Ljung", lag = 20)

    Box-Ljung test

data:  modelo.arma$residuals
X-squared = 2.1167, df = 20, p-value = 1

La hipotesis asociada a la prueba es la siguiente:

\[\begin{align*} H_0: \rho_1 = \rho_2 =0 \\ H_a: \textrm{Al menos un parámetro es distinto de 0} \end{align*}\]

Si p-valor > 0.05 NO se rechaza H0

Debido a que el p-valor (1) es mayor a 0.05, NO se rechaza H0, por lo tanto, NO hay estructura de autocorrelación en el modelo, lo cual lo hace un modelo válido.

Refinación del Modelo ARMA(1,1)

Mediante este proceso se busca observar si hay parámetros que no son significativos, lo cual significaría que no aportan al modelo en cuestión. Esto se realizará con la prueba “Coeftest”.

coeftest(modelo.arma, df = 67)

t test of coefficients:

          Estimate Std. Error t value  Pr(>|t|)    
ar1        0.75427    0.10985  6.8662 2.627e-09 ***
ma1       -0.19970    0.16191 -1.2334    0.2217    
intercept  8.61699    0.78141 11.0274 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#df = observaciones - parámetros
#df = 70 - 3

En esta prueba se observa que el MA es un parámetro NO significativo, por lo tanto, es necesario refinar el modelo. Esto se realizará asignando un cero a dicho parámetro NO significativo.

modelo.arma.ref <- arima(subocupacion, transform.pars = TRUE, order = c(1, 0, 1), fixed = c(NA,0,NA))
modelo.arma.ref

Call:
arima(x = subocupacion, order = c(1, 0, 1), transform.pars = TRUE, fixed = c(NA, 
    0, NA))

Coefficients:
         ar1  ma1  intercept
      0.6395    0     8.6172
s.e.  0.0895    0     0.6818

sigma^2 estimated as 4.442:  log likelihood = -151.78,  aic = 309.56

Correlograma de Modelo Refinado

par(mfrow = c(1,2))
acf(modelo.arma.ref$residuals, lag = 20)
pacf(modelo.arma.ref$residuals, lag = 20)

Debido a que se observa que todos las estructuras se encuentran dentro de las bandas de confianza, este será el nuevo modelo a utilizar. Es necesario mencionar que dicho modelo refinado se transforma en un modelo AR1, debido a que el MA no sale significativa.

Probar Estacionariedad

Finalmente, se realizará una prueba de estacionariedad mediante la obtención de las raíces del polinomio característico. En dicha prueba es necesario tomar en cuenta el valor de la raíz, debido a que si esta es mayor a 1, se considera que el modelo es estacionario.

pol.car.arma <-c(1,-modelo.arma.ref$coef[1])
raices <- polyroot(pol.car.arma)
abs(raices)

[1] 1.563666

Debido a que la raíz del polinomio característico es mayor a 1, se concluye que dicho modelo es estacionario.

2. Pronóstico de la Tasa de Subocupación

Con el objetivo de utilizar el mejor modelo adecuado para realizar un pronóstico óptimo, es necesario generar otro modelo en el cual se utilizará auto arima; una función automática que expone el mejor modelo a su consideración. Por lo tanto, dicho modelo estimado por auto arima, será comparado con el modelo ARMA (1,1) propuesto en la sección anterior.

Estimación Auto Arima

auto.sub <- auto.arima(subocupacion, d = 0, seasonal = F, ic = "bic")
auto.sub

Series: subocupacion 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1    mean
      0.6395  8.6172
s.e.  0.0895  0.6818

sigma^2 = 4.573:  log likelihood = -151.78
AIC=309.56   AICc=309.92   BIC=316.31

El modelo propuesto por la función auto arima arroja una sugerencia de modelo AR1.

Es necesario mencionar que en la sección anterior se encontró como modelo refinado un AR(1), asimismo, el modelo autoarima arrojó AR(1). Con el objetivo de realizar una comparación óptima entre modelos distintos, se decidió utilizar el modelo original (ARMA(1,1)) y no el refinado, con el fin de realizar una comparación con el modelo autoarima (AR(1)).

Posteriormente, se validara el modelo mediante un correlograma, en donde se identifique si dicho modelo sigue un proceso de Ruido Blanco, lo cual, significaría que no tiene autocorrelación, por lo tanto, esta sería estadísticamente cero.

Validación Modelo Auto Arima

par(mfrow = c(1,2))
acf(auto.sub$residuals, lag = 20)
pacf(auto.sub$residuals, lag = 20)

Se observa que el modelo AR(1) sigue un proceso de Ruido Blanco, por lo tanto, es un modelo adecuado.

En la generación de pronósticos se puede presentar el “error sobreajuste”, el cual quiere decir que los valores de los pronósticos es generado conforme a valores de “y”. Por lo tanto, se debe de dividir la muestra en grupos de entrenamiento y prueba. El conjunto de entrenamiento sirve para realizar un pronóstico que después se comparará con los valores de prueba.

Generación de grupos de entrenamiento y prueba

entrenamiento <- subset(sub.ts[,"sub"], end = length(sub.ts[,"sub"])- 14)
prueba <- subset(sub.ts[,"sub"], start = length(sub.ts[,"sub"])- 13)

Estimación de Modelos

arma11.entrenamiento <- Arima(entrenamiento, order = c(1,0,1))
ar1.entrenamiento <- Arima(entrenamiento, order = c(1,0,0))

Representación Gráfica

Pronóstico del modelo ARMA (1,1):

#ARMA(1,1)

arma11.entrenamiento %>%
  forecast(h=20) %>%
autoplot(xlab = "", ylab = "Tasa de subocupación (%)") + autolayer(prueba)

Pronóstico del modelo AR(1):

#AR(1)
ar1.entrenamiento %>%
  forecast(h=20) %>%
autoplot(xlab = "", ylab = "Tasa de subocupación (%)") + autolayer(prueba)

Con el objetivo de comparar cuál de los pronósticos es mejor, se realiza lo siguiente.

Comparación de ambos modelos mediante la prueba Diebold Mariano

La prueba Diebold Mariano, tiene como objetivo analizar si existe una diferencia sistemática entre los errores de los modelos, en dado caso de que sí, se considera mejor modelo el cual tenga un error menor.

sub.pron.arma11 <- forecast(arma11.entrenamiento, h = 20)
sub.pron.ar1 <- forecast(ar1.entrenamiento, h = 20)

Error del pronóstico

e.arma11 <- prueba - sub.pron.arma11$mean
e.ar1 <- prueba - sub.pron.ar1$mean

Error absoluto medio (MAE)

El Error Absoluto Medio se calcula con la siguiente fórmula:

\[\begin{align*} MAE = \frac{1}{p} \sum_{i=1}^{p}\left|y_{1}-f_{i} \right| \end{align*}\]

options(digits = 4)
MAE.arma11 <- sum(abs(e.arma11))/length(prueba)
MAE.arma11

[1] 4.14

MAE.ar1 <- sum(abs(e.ar1))/length(prueba)
MAE.ar1

[1] 4.059

En este caso los errores reflejan un valor alto debido a la gráfica observada, puesto que durante la pandemia se observaron resultados que no siguieron la tendencia de dicha tasa de subocupación, ya que este evento se relaciona con la teoría black swan, en donde dicho evento desató un gran impacto en el país.

Prueba Diebold Mariano

Debido a que no es suficiente observar el mejor modelo solo comparando los MAE, es necesario realizar una prueba que arroje si existe una diferencia sistemática entre los errores.

#Adaptación a rolling window (se debe hacer previo a la prueba DM)
farma11 <- function(x, h){forecast(Arima(x, order=c(1,0,1)), h=h)}
far1 <- function(x, h){forecast(Arima(x, order=c(1,0,0)), h=h)}
e.arma11.rw <- tsCV(sub.ts[,"sub"], farma11, h = 1, window = 56)
e.ar1.rw <- tsCV(sub.ts[,"sub"], far1, h = 1, window = 56)
dm.test(abs(e.arma11.rw), abs(e.ar1.rw),
        alternative = "two.sided", h = 1, power = 1)

    Diebold-Mariano Test

data:  abs(e.arma11.rw)abs(e.ar1.rw)
DM = 0.77, Forecast horizon = 1, Loss function power = 1, p-value = 0.4
alternative hypothesis: two.sided

La hipotesis asociada a la prueba es la siguiente: \[\begin{align*} H_0: \bar{d} = 0 \\ H_a: \bar{d} \neq 0 \end{align*}\]

Debido a que el p-valor > 0.05, NO se rechaza H0. Por lo tanto, no hay una diferencia sistemática entre los errores.

Sin embargo, tomando en cuenta que se tomó ARMA(1,1) solo con el objetivo de realizar una comparación entre lo propuesto por auto arima, y que con el modelo refinado se obtuvo un AR(1), se tomará como mejor modelo AR(1).

Pronóstico

pron.ar1 = forecast(auto.sub, h = 20)
pron.ar1

ar1.entrenamiento %>%
  forecast(h=20) %>%
autoplot(xlab = "", ylab = "Tasa de subocupación (%)") + autolayer(prueba)

CONCLUSIÓN

De acuerdo con el modelo AR(1), el pronóstico para el siguiente trimestre de la tasa de subocupación en México, es de 8.695%, se observa que este porcentaje va disminuyendo a lo largo del tiempo, hasta que transcurren aproximadamente 12 periodos (trimestres), que es cuando tiende a la media (8.617%).

REFERENCIAS BIBLIOGRÁFICAS

• Morales, F. (s.f.). Subempleo. https://economipedia.com/definiciones/subempleo.html
• Judisman, C. (s.f.). Conceptos y Definiciones en Relación con el Empleo, Desempleo y el Subempleo. https://estudiosdemograficosyurbanos.colmex.mx
  
• INEGI. (s.f.). Banco de información Económica (BIE). http://en.www.inegi.org.mx/app/indicadores/?tm=0&t=10000215#D10000215#D493912#D493911

ANEXOS

CÓDIGO RSTUDIO aquí