por:

Karla Isayuvi Amaro Estrada-1912032

Nevid Antonio Torres Lezama-1912046

Jose Miguel Manrique Velasco-1912040

Daniel Alfonso Hernández Gutiérrez-1912039

ANALSIS DEL SECTOR INDUSTRIAL

En México, durante el periodo comprendido de 2010 a 2020, la actividad industrial ha representado en promedio un 32% de la actividad económica total. Teniendo su participación más baja en el segundo trimestre del 2020 {ver gráfico 1}. La actividad industrial por si sola logra la absorción del 13% de la población ocupada, lo que representa, en términos absolutos, cerca de 12.5 millones de personas.

Fuente: Elaboración propia con base en BIE, INEGI.

La actividad industrial toatal puede dividirs en cuatro sub sectores: Minería; Generación, Transmisión y Distribución de Energía Eléctrica; Suministro de Agua y de Gas por Ductos al Consumidor Final; Construcción e Industrias Manufactureras. Esta última subdivisión, de Industrias Manufactureras, concentra el mayor dinamismo dentro de las actividades secundarias, por si sola, representa en promedio, el 52% de la actividad industrial total, siendo el 57% del valor de la actividad industrial total durante el 4to trimestre de 2020.

Fuente: Elaboración propia con base en BIE, INEGI.

Es por ello que variaciones que afecten la actividad manufacturera tendrán fuertes implicaciones en el comportamiento total de la actividad industrial. En este sentido, es relevante señalar la importancia del sector de los semiconductores dentro de la actividad manufacturera; las constantes innovaciones tecnológicas han tenido un fuerte impacto en la necesidad de las empresas, pues mientras hace unos años la principal necesidad de insumos recaía en metales brutos, hoy día el 80% de los costos de producción recae en desarrollos tecnológicos como los semiconductores, de tal suerte que la producción de los mismos podrá promover o afectar las cadenas de producción internacionales que emplean este importante insumo.

La industria automotriz mexicana es un claro ejemplo de la dinámica descrita anteriormente, pues la necesidad de insumos tecnológicos es esencial para su desarrollo. Además, una gran parte de la actividad económica de México se sostiene del crecimiento de la industria automotriz, por lo que, de manera indirecta, la actividad industrial depende de lo que suceda en el sector de los semiconductores, el cual ha sido afectado de manera importante por la pandemia.

En los primeros meses de la pandemia, la demanda de semiconductores cambio drásticamente; la industria automotriz dejó de ser el principal demandante de estos insumos, debido a la caída en las ventas de automóviles particulares y de transporte. Sin embargo, los proveedores de estos insumos lograron colocar su oferta en las industrias de televisores y computadoras, que aumentaron sus ventas dadas las nuevas dinámicas generadas por la pandemia. No obstante, a esto debemos sumar el contexto internacional, donde los principales proveedores han visto afectado su proceso de producción; por ejemplo, en EE.UU. las heladas sufridas durante la época invernal obligaron el cierre de empresas productoras de semiconductores, mientras que, en China, la producción se destinó principalmente al sector interno dadas las sanciones impuestas por EE.UU. en la guerra comercial promovida por el expresidente Donald Trump.

Fuente: Elaboración propia con base en datos de Edwin Tapia.

Así, una caída en la oferta de semiconductores a nivel mundial y la escaza oferta nacional de los mismos, combinado con una recuperación de las ventas del sector automotriz generará escasez de semiconductores para la industria automotriz, lo cual impactará negativamente la actividad industrial nacional. Esto lo podemos ver reflejado no solo en nuestro pronóstico de la actividad industrial: 98.22606, si no también en las noticias relacionadas al paro de actividades en diferentes plantas del sector automotriz.

REPORTE ESTADÍSTICO

Ventana de datos (Muestra): La ventana de tiempo utilizada para la estimación de nuestro modelo y posterior pronostico comprende de Ene-1993 a Feb-2021. Nuestra información es reportada de forma mensual.

Nuestras variables son: Actividad Industrial (AI), Producción de Vehículos Automóviles (PVA), Producción de Vehículos Camiones (PVC) e Indicador IMEG del Sector Manufacturero (IMEG). Es importante señalar que las bases de datos proporcionadas PVA y PVC están ajustados por estacionalidad.

ESTRUCTURA DE LOS DATOS

Modelo ARMAX Definición: El modelo ARMAX es un Proceso Autoregresivo de Media Móvil con Variables Exógenas. Esta forma de modelar es una variante del ARMA, donde estamos considerando que nuestro termino de error tiene su propia estructura de rezagos. En este caso tenemos que el modelo es multivariado,

Para la adecuada modelación se debe considerar que no hay estacionalidad en las series que utilizamos.

• Serie desestacionalizada de la PVA de México; periodicidad mensual, 1993-01 : 2021-02

• Serie desestacionalizada de la PVC de México; periodicidad mensual, 1993-01 : 2021-02

• Serie desestacionalizada de la AI de México; periodicidad mensual, 1993-01 : 2021-02

• Serie desestacionalizada de la AI de México; periodicidad mensual, 1993-01 : 2021-01

• Serie desestacionalizada IMEG; periodicidad mensual, 2005-01 : 2021-02

Paso 1: Elegimos la variable independiente que mejor explique el crecimiento de la actividad industrial de México (Y)

i<- NULL
lag_Xs <- c()

for(i in 1:12){
  if (i==1){   
    lag_Xs<-cbind(X1,lag.xts(X1,k=i))
  } else {
    lag_Xs <-cbind(lag_Xs,lag.xts(X1,k=i))
  }
}
for(i in 1:12){
  if (i==1){   
    lag_Xs<-cbind(lag_Xs,X2,lag.xts(X2,k=i))
  } else {
    lag_Xs <-cbind(lag_Xs,lag.xts(X2,k=i))
  }
}


colnames(lag_Xs)[1:13]<- paste("prod_auto", 0:12, sep = "") 
colnames(lag_Xs)[14:26]<- paste("prod_camio", 0:12, sep = "")
act_ind_Xs<-data.frame(Y,lag_Xs)

ggcorrplot(cor(act_ind_Xs[13:NROW(act_ind_Xs),]),lab=T, lab_size = 2 , type = "upper")

Para este primer análisis se empleó la base de producción de automóviles y camiones en México, pues como mencionamos en el análisis de la industria, es el sector más dinámico dentro de la actividad económica industrial. Podemos ver en nuestro mapa de correlación que las variables explicativas son contemporáneas, es decir emplearemos dentro de nuestras variables explicativas (X) la actividad en el periodo 0 (sin rezagos).

Por ello tenemos que nuestro modelo tendrá la forma:

\(Y_t=ϕ(L) X_t+θ(L) ε_t\)

.

Paso 2: Estimamos la regresión entre X & Y

xreg <- lag_Xs[, c(1,14)] # Tomamos la variables contemporáneas

dummy <- ifelse( mean(Y)-2*sd(Y) < Y & Y < mean(Y)-sd(Y), 0.8,
              ifelse(Y <= (mean(Y)-2*sd(Y)), 1,0))
dum_date <- data.frame(Y, dummy , 
                       dates =  seq(as.Date("1993-02-01"),
                                    as.Date("2021-01-01"),
                                    by = "month"))

xreg1 <- cbind(xreg,dummy )


m1 <- lm(Y ~ xreg1)
summary(m1)

Call:
lm(formula = Y ~ xreg1)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.9275  -0.7061  -0.0891   0.6035   7.6156 

Coefficients:
                  Estimate Std. Error t value             Pr(>|t|)    
(Intercept)       0.222867   0.079617   2.799              0.00542 ** 
xreg1prod_auto0   0.011635   0.002888   4.029            0.0000696 ***
xreg1prod_camio0  0.045090   0.004183  10.780 < 0.0000000000000002 ***
xreg1            -5.987405   0.553487 -10.818 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.433 on 332 degrees of freedom
Multiple R-squared:  0.6159,    Adjusted R-squared:  0.6124 
F-statistic: 177.4 on 3 and 332 DF,  p-value: < 0.00000000000000022

fit<-as.matrix(fitted(m1))
plot.ts(coredata(Y),col="red") 
lines(fit,col="blue")

library(urca)

residuales_m1=m1$residuals
residualPlot(m1)

DF_m1 =ur.df(residuales_m1)

DF_m1@teststat

               tau1
statistic -13.19334

Paso 3: Revisamos la función de autocorrelación de los residuos para encontrar un modelo ARMA que explique los residuos

acf2(residuals(m1))

      [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  -0.14 0.06 -0.04 0.03 0.05 -0.08 -0.03 0.02 -0.03  0.02  0.04 -0.07  0.09
PACF -0.14 0.04 -0.03 0.01 0.06 -0.07 -0.06 0.02 -0.03  0.00  0.06 -0.07  0.07
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.07  0.07  0.01  0.00  0.04  0.04  0.04  0.01 -0.05  0.03 -0.07  0.06
PACF -0.03  0.03  0.04  0.01  0.03  0.06  0.05  0.02 -0.04  0.02 -0.07  0.06
     [,26] [,27] [,28] [,29]
ACF  -0.01  0.02 -0.05  0.02
PACF  0.00  0.03 -0.06  0.01

Con nuestro plot de residuos podemos considerar un MA(24)

Paso 4: Estimamos nuestro modelo ARMAX

y_am1 <- sarima(Y,1,0,0,
                 fixed = c(NA,
                           NA,
                           rep(NA,3))
                 ,xreg=xreg1)

initial  value 0.353749 
iter   2 value 0.343527
iter   3 value 0.341720
iter   4 value 0.338104
iter   5 value 0.338081
iter   6 value 0.338060
iter   7 value 0.338059
iter   8 value 0.338059
iter   9 value 0.338059
iter   9 value 0.338059
iter   9 value 0.338059
final  value 0.338059 
converged
initial  value 0.337896 
iter   2 value 0.337894
iter   3 value 0.337893
iter   4 value 0.337893
iter   5 value 0.337893
iter   5 value 0.337893
iter   5 value 0.337893
final  value 0.337893 
converged

# y_am11 <- sarima(Y,1,0,1,
#                  fixed = c(NA,
#                            NA,
#                            NA,
#                            rep(NA,3))
#                  ,xreg=xreg1)

# y_am1         # 3.574219 sin media / 3.549376 con media
# y_am11        # 3.574115 sin media / 3.551223 con media   # MA no significativo

acf2(y_am1[["fit"]][["residuals"]])

     [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  0.01 0.07 -0.01 0.02 0.05 -0.07 -0.05 -0.01 -0.04  0.02  0.04 -0.05  0.07
PACF 0.01 0.07 -0.01 0.02 0.05 -0.08 -0.05  0.01 -0.04  0.02  0.05 -0.06  0.06
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.05  0.07  0.03     0  0.04  0.06  0.06  0.01 -0.04  0.01 -0.06  0.04
PACF -0.04  0.05  0.03     0  0.02  0.07  0.05  0.00 -0.03  0.01 -0.06  0.06
     [,26] [,27] [,28] [,29]
ACF      0  0.01 -0.05  0.01
PACF     0  0.02 -0.05  0.02

mean(y_am1[["fit"]][["residuals"]])

[1] 0.000787159

Paso 5: Forecast. Pronosticamos para cada una de nuestras variables

xpred1 <- cbind(X1c[338], X2c[338], 0 )

Y1_Pred <- sarima.for(Y,1,0,0,
                         fixed = c(NA,
                                   NA,
                                   rep(NA,3)),
                         xreg=xreg1,
                         newxreg = xpred1,
                         n.ahead = 1)

ARMAX_Feb21 <- as.numeric(Y1_Pred$pred)
data$act_ind[337]*exp(ARMAX_Feb21/100)

[1] 96.98933

Paso 1: Elegimos la variable independiente que mejor explique el crecimienmto de la actividad industrial de México (Y)

i<- NULL
lag_2Xs <- c()

for(i in 1:12){
  if (i==1){   
    lag_2Xs<-cbind(Z1,lag.xts(Z1,k=i))
  } else {
    lag_2Xs <-cbind(lag_2Xs,lag.xts(Z1,k=i))
  }
}
for(i in 1:12){
  if (i==1){   
    lag_2Xs<-cbind(lag_2Xs,Z2,lag.xts(Z2,k=i))
  } else {
    lag_2Xs <-cbind(lag_2Xs,lag.xts(Z2,k=i))
  }
}


colnames(lag_2Xs)[1:13]<- paste("pmi", 0:12, sep = "") 
colnames(lag_2Xs)[14:26]<- paste("conf_man", 0:12, sep = "")
act_ind_2Xs<-data.frame(Y,lag_2Xs)

ggcorrplot(cor(act_ind_2Xs[13:NROW(act_ind_2Xs),]),lab=T, lab_size = 2 , type = "upper")

Se observa que las variables explicativas son contemporáneas

Paso 2: Estimamos la regresion entre X & Y

zreg <- lag_2Xs[, c(1,14)] # Tomamos la variables contemporáneas

dummy2 <- ifelse( mean(Y)-2*sd(Y) < Y & Y < mean(Y)-sd(Y), 0.8,
              ifelse(Y <= (mean(Y)-2*sd(Y)), 1,0))
sum(dummy2 != 0)/NROW(dummy2)*100

[1] 1.5625

dum2_date <- data.frame(Y, dummy2 , 
                       dates =  seq(as.Date("2005-02-01"),
                                    as.Date("2021-01-01"),
                                    by = "month"))

zreg1 <- cbind(zreg,dummy2 )


m2 <- lm(Y ~ zreg1)
summary(m2)

Call:
lm(formula = Y ~ zreg1)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.9009  -0.7245  -0.1369   0.5204  12.9664 

Coefficients:
                Estimate Std. Error t value             Pr(>|t|)    
(Intercept)      0.23132    0.14234   1.625             0.105814    
zreg1pmi0        0.16182    0.04502   3.595             0.000415 ***
zreg1conf_man0   0.19330    0.05567   3.472             0.000640 ***
zreg1act_ind   -13.10040    1.40991  -9.292 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.957 on 188 degrees of freedom
Multiple R-squared:  0.4955,    Adjusted R-squared:  0.4875 
F-statistic: 61.56 on 3 and 188 DF,  p-value: < 0.00000000000000022

fit<-as.matrix(fitted(m2))
plot.ts(Y,col="red") 
lines(fit,col="blue")

#Realizamos la prueba de raíz unitaria DFuller

residuales_m2=m2$residuals
residualPlot(m2)

# Con la gráfica podemos intuir que los datos no son estacionarios

DF_m2 =ur.df(residuales_m2)
summary(DF_m2)

############################################### 
# 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 
-11.4558  -0.8291  -0.1885   0.5912  10.7222 

Coefficients:
           Estimate Std. Error t value            Pr(>|t|)    
z.lag.1    -1.33073    0.10737 -12.393 <0.0000000000000002 ***
z.diff.lag  0.18566    0.07165   2.591              0.0103 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.908 on 188 degrees of freedom
Multiple R-squared:  0.5763,    Adjusted R-squared:  0.5718 
F-statistic: 127.8 on 2 and 188 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -12.3935 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

DF_m2@teststat

               tau1
statistic -12.39346

Se corrobora que los datos son estacionarios, pues nuestro T estadistico en la prueba de DFuller es mayor en valor absoluto a 10

Paso 3: Revisamos la función de utocorrelación de los residuos para encontrar un modelo ARMA que explique los residuos

acf2(residuals(m2))

      [,1]  [,2] [,3] [,4] [,5]  [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
ACF  -0.12 -0.17 0.09 0.05 0.00 -0.01 0.01 -0.01 0.01 -0.08  0.05  0.00  0.04
PACF -0.12 -0.19 0.04 0.04 0.04  0.01 0.01 -0.02 0.00 -0.09  0.03 -0.02  0.06
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
ACF  -0.01 -0.05  0.00 -0.01 -0.01     0 -0.02  0.01 -0.02  0.01 -0.09
PACF  0.00 -0.03 -0.03 -0.04 -0.03     0 -0.03  0.02 -0.03  0.02 -0.11

Con nuestro plot de residuos podemos considerar un MA(2)

Paso 4: Estimamos nuestro modelo ARMAX

y2_am2 <- sarima(Y,0,0,2,
                 xreg=zreg1)

initial  value 0.660806 
iter   2 value 0.640730
iter   3 value 0.635339
iter   4 value 0.633199
iter   5 value 0.633097
iter   6 value 0.633084
iter   7 value 0.633084
iter   8 value 0.633084
iter   8 value 0.633084
iter   8 value 0.633084
final  value 0.633084 
converged
initial  value 0.633390 
iter   2 value 0.633388
iter   3 value 0.633386
iter   4 value 0.633386
iter   5 value 0.633386
iter   5 value 0.633386
iter   5 value 0.633386
final  value 0.633386 
converged

acf2(y_am1[["fit"]][["residuals"]])

     [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  0.01 0.07 -0.01 0.02 0.05 -0.07 -0.05 -0.01 -0.04  0.02  0.04 -0.05  0.07
PACF 0.01 0.07 -0.01 0.02 0.05 -0.08 -0.05  0.01 -0.04  0.02  0.05 -0.06  0.06
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.05  0.07  0.03     0  0.04  0.06  0.06  0.01 -0.04  0.01 -0.06  0.04
PACF -0.04  0.05  0.03     0  0.02  0.07  0.05  0.00 -0.03  0.01 -0.06  0.06
     [,26] [,27] [,28] [,29]
ACF      0  0.01 -0.05  0.01
PACF     0  0.02 -0.05  0.02

mean(y_am1[["fit"]][["residuals"]])

[1] 0.000787159

Paso 5: Forecast. Pronosticamos para cada una de nuetras variables

zpred1 <- cbind(Z1c[194], Z2c[194], 0 )

Y2_Pred <- sarima.for(Y,0,0,2,
                         xreg= zreg1,
                         newxreg = zpred1,
                         n.ahead = 1)

ARMAX2_Feb21 <- as.numeric(Y2_Pred$pred)
data2$act_ind[193]*exp(ARMAX2_Feb21/100)

[1] 98.22606

Paso 1: Elegimos la variable independiente que mejor explique el crecimienmto de la actividad industrial de México (Y)

i<- NULL
lag_3Xs <- c()

for(i in 1:12){
  if (i==1){   
    lag_3Xs<-cbind(B1,lag.xts(B1,k=i))
  } else {
    lag_3Xs <-cbind(lag_3Xs,lag.xts(B1,k=i))
  }
}
for(i in 1:12){
  if (i==1){   
    lag_3Xs<-cbind(lag_3Xs,B2,lag.xts(B2,k=i))
  } else {
    lag_3Xs <-cbind(lag_3Xs,lag.xts(B2,k=i))
  }
}


colnames(lag_3Xs)[1:13]<- paste("pmi", 0:12, sep = "") 
colnames(lag_3Xs)[14:26]<- paste("ind_pedman", 0:12, sep = "")
act_ind_3Xs<-data.frame(Y,lag_3Xs)

ggcorrplot(cor(act_ind_3Xs[13:NROW(act_ind_3Xs),]),lab=T, lab_size = 2 , type = "upper")

Se observa que las variables explicativas son contemporáneas

Paso 2: Estimamos la regresion entre X & Y

breg <- lag_3Xs[, c(1,14)] # Tomamos la variables contemporáneas

# dummy2 <- ifelse( mean(Y)-2*sd(Y) < Y & Y < mean(Y)-sd(Y), 0.8,
#               ifelse(Y <= (mean(Y)-2*sd(Y)), 1,0))
# sum(dummy2 != 0)/NROW(dummy2)
# dum2_date <- data.frame(Y, dummy2 , 
#                        dates =  seq(as.Date("2005-02-01"),
#                                     as.Date("2021-01-01"),
#                                     by = "month"))

breg1 <- cbind(breg,dummy2)


m3 <- lm(Y ~ breg1)
summary(m3)

Call:
lm(formula = Y ~ breg1)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.5168  -0.8303  -0.1240   0.5430   9.0658 

Coefficients:
                  Estimate Std. Error t value             Pr(>|t|)    
(Intercept)        0.20160    0.13344   1.511              0.13251    
breg1pmi0          0.11909    0.04278   2.784              0.00592 ** 
breg1ind_pedman0   0.33009    0.05223   6.320        0.00000000185 ***
breg1act_ind     -12.74464    1.28979  -9.881 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.833 on 188 degrees of freedom
Multiple R-squared:  0.5572,    Adjusted R-squared:  0.5502 
F-statistic: 78.87 on 3 and 188 DF,  p-value: < 0.00000000000000022

fit<-as.matrix(fitted(m3))
plot.ts(coredata(Y),col="red")
lines(fit,col="blue")

#Realizamos la prueba de raíz unitaria DFuller

residuales_m3=m3$residuals
residualPlot(m3)

# Con la gráfica podemos intuir que los datos no son estacionarios

DF_m3 =ur.df(residuales_m3)
summary(DF_m3)

############################################### 
# 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 
-9.6814 -0.7910 -0.1513  0.5523  8.8485 

Coefficients:
           Estimate Std. Error t value             Pr(>|t|)    
z.lag.1    -1.52977    0.11359 -13.467 < 0.0000000000000002 ***
z.diff.lag  0.20612    0.07128   2.892              0.00428 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.726 on 188 degrees of freedom
Multiple R-squared:  0.6497,    Adjusted R-squared:  0.646 
F-statistic: 174.3 on 2 and 188 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -13.4674 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

DF_m3@teststat

               tau1
statistic -13.46745

Paso 3: Revisamos la función de autocorrelación de los residuos para encontrar un modelo ARMA que explique los residuos

acf2(residuals(m3))

      [,1]  [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
ACF  -0.27 -0.12  0.07 0.06 -0.04 0.04 -0.03 -0.01 0.03 -0.08  0.05 -0.06  0.05
PACF -0.27 -0.21 -0.03 0.06  0.01 0.06 -0.01 -0.01 0.02 -0.08  0.01 -0.08  0.03
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
ACF   0.01 -0.01  0.00 -0.01 -0.01  0.02 -0.02 -0.05  0.05  0.01 -0.10
PACF  0.03  0.01  0.02 -0.02 -0.02  0.00 -0.03 -0.05  0.00  0.03 -0.08

Con nuestro plot de residuos podemos considerar un MA(1)

Paso 4: Estimamos nuestro modelo ARMAX

y3_am1 <- sarima(Y,0,0,1,
                 xreg=breg1)

initial  value 0.595555 
iter   2 value 0.547202
iter   3 value 0.536239
iter   4 value 0.531976
iter   5 value 0.531546
iter   6 value 0.531524
iter   7 value 0.531521
iter   8 value 0.531521
iter   9 value 0.531521
iter   9 value 0.531521
iter   9 value 0.531521
final  value 0.531521 
converged
initial  value 0.531694 
iter   2 value 0.531693
iter   3 value 0.531693
iter   4 value 0.531693
iter   4 value 0.531693
iter   4 value 0.531693
final  value 0.531693 
converged

acf2(y_am1[["fit"]][["residuals"]])

     [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  0.01 0.07 -0.01 0.02 0.05 -0.07 -0.05 -0.01 -0.04  0.02  0.04 -0.05  0.07
PACF 0.01 0.07 -0.01 0.02 0.05 -0.08 -0.05  0.01 -0.04  0.02  0.05 -0.06  0.06
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.05  0.07  0.03     0  0.04  0.06  0.06  0.01 -0.04  0.01 -0.06  0.04
PACF -0.04  0.05  0.03     0  0.02  0.07  0.05  0.00 -0.03  0.01 -0.06  0.06
     [,26] [,27] [,28] [,29]
ACF      0  0.01 -0.05  0.01
PACF     0  0.02 -0.05  0.02

mean(y_am1[["fit"]][["residuals"]])

[1] 0.000787159

Paso 5: Forecast. Pronosticamos para cada una de nuetras variables

bpred1 <- cbind(B1c[194], B2c[194], 0 )

Y3_Pred <- sarima.for(Y,0,0,1,
                         xreg= breg1,
                         newxreg = bpred1,
                         n.ahead = 1)

ARMAX3_Feb21 <- as.numeric(Y3_Pred$pred)
data2$act_ind[193]*exp(ARMAX3_Feb21/100)

[1] 97.2746

Paso 1: Elegimos la variable independiente que mejor explique el crecimienmto de la actividad industrial de México (Y)

i<- NULL
lag_4Xs <- c()

for(i in 1:12){
  if (i==1){   
    lag_4Xs<-cbind(D1,lag.xts(D1,k=i))
  } else {
    lag_4Xs <-cbind(lag_4Xs,lag.xts(D1,k=i))
  }
}
for(i in 1:12){
  if (i==1){   
    lag_4Xs<-cbind(lag_4Xs,D2,lag.xts(D2,k=i))
  } else {
    lag_4Xs <-cbind(lag_4Xs,lag.xts(D2,k=i))
  }
}


colnames(lag_4Xs)[1:13]<- paste("pmi", 0:12, sep = "") 
colnames(lag_4Xs)[14:26]<- paste("imp", 0:12, sep = "")
act_ind_4Xs<-data.frame(Y,lag_4Xs)

ggcorrplot(cor(act_ind_4Xs[13:NROW(act_ind_4Xs),]),lab=T, lab_size = 2 , type = "upper")

Se observa que las variables explicativas son contemporáneas

Paso 2: Estimamos la regresion entre X & Y

dreg <- lag_4Xs[, c(1,14)] # Tomamos la variables contemporáneas

# dummy2 <- ifelse( mean(Y)-2*sd(Y) < Y & Y < mean(Y)-sd(Y), 0.8,
#               ifelse(Y <= (mean(Y)-2*sd(Y)), 1,0))
# sum(dummy2 != 0)/NROW(dummy2)
# dum2_date <- data.frame(Y, dummy2 , 
#                        dates =  seq(as.Date("2005-02-01"),
#                                     as.Date("2021-01-01"),
#                                     by = "month"))

dreg1 <- cbind(dreg,dummy2)


m4 <- lm(Y ~ dreg1)
summary(m4)

Call:
lm(formula = Y ~ dreg1)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5306  -0.8176  -0.0243   0.5129  12.9501 

Coefficients:
              Estimate Std. Error t value             Pr(>|t|)    
(Intercept)    0.20013    0.14350   1.395             0.164766    
dreg1pmi0      0.17374    0.04465   3.892             0.000138 ***
dreg1imp0      0.05184    0.01651   3.139             0.001968 ** 
dreg1act_ind -14.01036    1.36123 -10.292 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.968 on 188 degrees of freedom
Multiple R-squared:  0.4899,    Adjusted R-squared:  0.4818 
F-statistic: 60.19 on 3 and 188 DF,  p-value: < 0.00000000000000022

fit<-as.matrix(fitted(m4))
plot.ts(Y,col="red")
lines(fit,col="blue")

#Realizamos la prueba de raíz unitaria DFuller

residuales_m4=m4$residuals
residualPlot(m4)

# Con la gráfica podemos intuir que los datos no son estacionarios

DF_m4 =ur.df(residuales_m4)
summary(DF_m4)

############################################### 
# 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 
-11.5904  -0.7763  -0.0916   0.5182  10.3774 

Coefficients:
           Estimate Std. Error t value             Pr(>|t|)    
z.lag.1    -1.42251    0.10918 -13.029 < 0.0000000000000002 ***
z.diff.lag  0.21152    0.07124   2.969              0.00338 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.893 on 188 degrees of freedom
Multiple R-squared:  0.6056,    Adjusted R-squared:  0.6014 
F-statistic: 144.3 on 2 and 188 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -13.0291 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

DF_m4@teststat

               tau1
statistic -13.02908

Paso 3: Revisamos la función de autocorrelación de los residuos para encontrar un modelo ARMA que explique los residuos

acf2(residuals(m4))

      [,1]  [,2] [,3]  [,4] [,5] [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
ACF  -0.17 -0.17 0.10  0.00 0.01 0.02 0.00 -0.04 0.03 -0.08  0.04  0.02  0.00
PACF -0.17 -0.21 0.03 -0.01 0.04 0.03 0.03 -0.03 0.02 -0.10  0.02 -0.01  0.03
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
ACF   0.04 -0.05 -0.02  0.00  0.01  0.02 -0.03 -0.02 -0.03  0.06 -0.06
PACF  0.05 -0.02 -0.02 -0.03 -0.01  0.03 -0.03 -0.02 -0.06  0.04 -0.05

Con nuestro plot de residuos podemos considerar un MA(2)

Paso 4: Estimamos nuestro modelo ARMAX

y4_am2 <- sarima(Y,0,0,2,
                 xreg=dreg1)

initial  value 0.666340 
iter   2 value 0.634152
iter   3 value 0.629615
iter   4 value 0.628139
iter   5 value 0.628058
iter   6 value 0.628053
iter   7 value 0.628053
iter   8 value 0.628053
iter   8 value 0.628053
iter   8 value 0.628053
final  value 0.628053 
converged
initial  value 0.628397 
iter   2 value 0.628395
iter   3 value 0.628393
iter   4 value 0.628392
iter   5 value 0.628392
iter   5 value 0.628392
iter   5 value 0.628392
final  value 0.628392 
converged

acf2(y_am1[["fit"]][["residuals"]])

     [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  0.01 0.07 -0.01 0.02 0.05 -0.07 -0.05 -0.01 -0.04  0.02  0.04 -0.05  0.07
PACF 0.01 0.07 -0.01 0.02 0.05 -0.08 -0.05  0.01 -0.04  0.02  0.05 -0.06  0.06
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.05  0.07  0.03     0  0.04  0.06  0.06  0.01 -0.04  0.01 -0.06  0.04
PACF -0.04  0.05  0.03     0  0.02  0.07  0.05  0.00 -0.03  0.01 -0.06  0.06
     [,26] [,27] [,28] [,29]
ACF      0  0.01 -0.05  0.01
PACF     0  0.02 -0.05  0.02

mean(y_am1[["fit"]][["residuals"]])

[1] 0.000787159

Paso 5: Forecast. Pronosticamos para cada una de nuetras variables

dpred1 <- cbind(D1c[194], D2c[194], 0 )

Y4_Pred <- sarima.for(Y,0,0,2,
                         xreg= dreg1,
                         newxreg = dpred1,
                         n.ahead = 1)

ARMAX4_Feb21 <- as.numeric(Y4_Pred$pred)
data2$act_ind[193]*exp(ARMAX4_Feb21/100)

[1] 97.73581