Modelo ARIMAX
Roberto Ryder López
10/04/2023
Los modelos ARIMAX
library(readxl)
library(tseries)
library(stargazer)
library(lmtest)
library(car)
library(strucchange)
library(TSstudio)
library(forecast)# Base de datos
data_arimax_mex <- read_excel("data_arimax_mex.xlsx",
sheet = "data4")
head(data_arimax_mex)Se declara la serie de tiempo.
# Tasa de crecimiento logaritmica del consumo privado
tlcp <- ts(data_arimax_mex[ , 2],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
ts_plot(tlcp,
title = "Consumo privado de México 1994.1-2022.4",
Ytitle = "Tasa de crecimiento logaritmica",
Xtitle = "Fuente: INEGI",
slider = TRUE)# Tasa de crecimiento logaritmica del PIB
tly <- ts(data_arimax_mex[ , 3],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
plot(tly)
# Tasa de crecimiento logaritmica de la inversión privada
tli <- ts(data_arimax_mex[ , 4],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
plot(tli)
# Tasa de crecimiento logaritmica de la inflacion
tlipc <- ts(data_arimax_mex[ , 5],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
plot(tlipc)
# Dummy 1
du1 <- ts(data_arimax_mex[ , 6],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
# Dummy 2
du2 <- ts(data_arimax_mex[ , 7],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
# Dummy 3
du3 <- ts(data_arimax_mex[ , 8],
start = c(1994, 1),
end = c(2023, 4),
freq = 4)
# Descomposicion
ts_decompose(tlcp)1 Estacionariedad
Para confirmar la estacionaridad de las series se usan las pruebas de Dickey-Fuller (ADF) y Phillips-Perron (PP) para raíz unitaria, donde:
- \(H_0\): La serie no es estacionaria. \(p.value \geq 0.05\)
- \(H_1\): La serie es estacionaria. \(p.value \leq 0.05\)
# Prueba de Dickey-Fuller (ADF)
tlcpn <- na.remove(tlcp)
adf.test(tlcpn)##
## Augmented Dickey-Fuller Test
##
## data: tlcpn
## Dickey-Fuller = -3.6959, Lag order = 4, p-value = 0.02762
## alternative hypothesis: stationarytlyn <- na.remove(tly)
adf.test(tlyn)##
## Augmented Dickey-Fuller Test
##
## data: tlyn
## Dickey-Fuller = -3.728, Lag order = 4, p-value = 0.02485
## alternative hypothesis: stationarytlin <- na.remove(tli)
adf.test(tlin)##
## Augmented Dickey-Fuller Test
##
## data: tlin
## Dickey-Fuller = -3.5291, Lag order = 4, p-value = 0.04265
## alternative hypothesis: stationarytlipcn <- na.remove(tlipc)
adf.test(tlipcn)##
## Augmented Dickey-Fuller Test
##
## data: tlipcn
## Dickey-Fuller = -2.0379, Lag order = 4, p-value = 0.5608
## alternative hypothesis: stationary # Primeras diferencias
dtlipcn <- diff(tlipcn)
adf.test(dtlipcn)##
## Augmented Dickey-Fuller Test
##
## data: dtlipcn
## Dickey-Fuller = -10.739, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary# Prueba de Phillips-Perron (PP)
pp.test(tlcpn)##
## Phillips-Perron Unit Root Test
##
## data: tlcpn
## Dickey-Fuller Z(alpha) = -43.163, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarypp.test(tlyn)##
## Phillips-Perron Unit Root Test
##
## data: tlyn
## Dickey-Fuller Z(alpha) = -45.98, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarypp.test(tlin)##
## Phillips-Perron Unit Root Test
##
## data: tlin
## Dickey-Fuller Z(alpha) = -34.401, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarypp.test(dtlipcn)##
## Phillips-Perron Unit Root Test
##
## data: dtlipcn
## Dickey-Fuller Z(alpha) = -48.089, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryTodas las series son estacionarias, en el caso de la inflación fue necesario aplicar la primera diferencia.
2 Diseño muestral
A continuación se crea la muestra para la estimación del modelo ARIMAX y la muestra para la proyección. Recuerda que la muestra exogena se divide en dos partes con este proposito.
# Endogena
estacionaria <- cbind(tlcp)
muestra_end <- window(estacionaria,
start = c(1994, 1),
end = c(2022,4),
freq = 4)
# Exogena
muestra_exo <- ts.intersect(tly,
tli,
tlipcn,
du1,
du2,
du3)
# Exogena para la estimacion
muestra_exo1 <- window(muestra_exo,
start = c(1994, 1),
end = c(2022, 4),
freq = 4)
# Exogena para la proyeccion
muestra_exo2 <- window(muestra_exo,
start = c(2023, 1),
end = c(2023, 4),
freq = 4)3 Correlograma
Se identifican los componentes p (rezagos) y q (medias móviles) con ayuda del correlograma.
# Correlograma
ts_cor(tlcpn, lag = 13)4 Estimación
arimax <- arima(muestra_end,
order = c(1, 0, 2),
xreg = muestra_exo1,
method = "ML")
summary(arimax)##
## Call:
## arima(x = muestra_end, order = c(1, 0, 2), xreg = muestra_exo1, method = "ML")
##
## Coefficients:
## ar1 ma1 ma2 intercept tly tli tlipcn du1
## 0.5335 0.0634 1.0000 -0.0085 0.9326 0.0175 0.0175 -0.0386
## s.e. 0.1104 0.0225 0.0622 0.0345 0.0494 0.0196 0.0094 0.0530
## du2 du3
## 0.0608 -0.0445
## s.e. 0.0526 0.0573
##
## sigma^2 estimated as 0.003887: log likelihood = 152.8, aic = -283.59
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0006078931 0.0623444 0.04807442 30.60329 52.79909 0.4077034
## ACF1
## Training set 0.1508565coeftest(arimax)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.5334815 0.1103525 4.8343 1.336e-06 ***
## ma1 0.0633842 0.0225391 2.8122 0.004921 **
## ma2 0.9999809 0.0621954 16.0781 < 2.2e-16 ***
## intercept -0.0085386 0.0344500 -0.2479 0.804246
## tly 0.9326386 0.0494351 18.8659 < 2.2e-16 ***
## tli 0.0174513 0.0195707 0.8917 0.372550
## tlipcn 0.0174631 0.0094230 1.8532 0.063847 .
## du1 -0.0385835 0.0530320 -0.7276 0.466889
## du2 0.0607508 0.0526106 1.1547 0.248202
## du3 -0.0444822 0.0572632 -0.7768 0.437275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15 Estimación 2
# Endogena
estacionaria <- cbind(tlcp)
muestra_end <- window(estacionaria,
start = c(1994, 1),
end = c(2022,4),
freq = 4)
# Exogena
muestra_exo <- ts.intersect(tly,
tlipcn)
# Exogena para la estimacion
muestra_exo1 <- window(muestra_exo,
start = c(1994, 1),
end = c(2022, 4),
freq = 4)
# Exogena para la proyeccion
muestra_exo2 <- window(muestra_exo,
start = c(2023, 1),
end = c(2023, 4),
freq = 4)
# Estimacion del modelo
arimax <- arima(muestra_end,
order = c(1, 0, 2),
xreg = muestra_exo1,
method = "ML")
summary(arimax)##
## Call:
## arima(x = muestra_end, order = c(1, 0, 2), xreg = muestra_exo1, method = "ML")
##
## Coefficients:
## ar1 ma1 ma2 intercept tly tlipcn
## 0.5221 0.0734 1.0000 -0.0110 0.9719 0.0171
## s.e. 0.1042 0.0228 0.0493 0.0328 0.0291 0.0089
##
## sigma^2 estimated as 0.004001: log likelihood = 151.13, aic = -288.25
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0006228784 0.0632556 0.05036652 35.00961 57.67578 0.4271419
## ACF1
## Training set 0.1546356coeftest(arimax)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.5220899 0.1042373 5.0087 5.481e-07 ***
## ma1 0.0734152 0.0227604 3.2256 0.001257 **
## ma2 0.9999939 0.0492756 20.2939 < 2.2e-16 ***
## intercept -0.0110108 0.0327574 -0.3361 0.736772
## tly 0.9718953 0.0291434 33.3488 < 2.2e-16 ***
## tlipcn 0.0170575 0.0089403 1.9079 0.056400 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1stargazer(arimax, type = "text")##
## =============================================
## Dependent variable:
## ---------------------------
## muestra_end
## ---------------------------------------------
## ar1 0.522***
## (0.104)
##
## ma1 0.073***
## (0.023)
##
## ma2 1.000***
## (0.049)
##
## intercept -0.011
## (0.033)
##
## tly 0.972***
## (0.029)
##
## tlipcn 0.017*
## (0.009)
##
## ---------------------------------------------
## Observations 116
## Log Likelihood 151.126
## sigma2 0.004
## Akaike Inf. Crit. -288.251
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01El modelo final es:
\[\widehat{tlcp}_t = - \widehat{\beta}_0 + ~ AR_1 ~ tlcp_{t-1} + ~ MA_1 + ~ MA_2 ~ + \beta_1 ~ tly + \beta_2 ~ tlipcn \] \[\widehat{tlcp}_t = - ~\widehat{0.01} ~ + \widehat{0.52} ~ tlcp_{t-1} + \widehat{0.07} ~ \epsilon_{t-1} + \widehat{1.00} ~ \epsilon_{t-2} ~ + ~ 0.97~ tly + ~ 0.02 ~ tlipcn\]
6 Evaluación
6.1 Normalidad
Es necesario que los residuos se comporten como ruido blanco, para ellos usamos la prueba de Dickey-Fuller aumentada, donde:
- \(H_0\): Los residuos no se comportan como ruido blanco. \(p.value \geq 0.05\)
- \(H_1\): Los residuos se comportan como ruido blanco. \(p.value \leq 0.05\)
# Prueba ADF para
residuals <- resid(arimax)
adf.test(residuals)##
## Augmented Dickey-Fuller Test
##
## data: residuals
## Dickey-Fuller = -4.779, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary6.2 Autocorrelación
Para determinar que un grupo cualquiera de autocorrelaciones sea distinta de cero, es decir, para comprobar si una serie de observaciones en un período de tiempo específico son aleatorias e independientes, se usa la prueba de Ljung-Box, donde:
- \(H_0\): No hay autocorrelacion, los residuos se distribuyen de forma independiente. \(p.value \geq 0.05\)
- \(H_1\): Si hay autocorrelacion, los residuos no se distribuyen de forma independiente \(p.value \leq 0.05\)
# Prueba Ljung-Box
checkresiduals(arimax)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2) with non-zero mean
## Q* = 17.775, df = 5, p-value = 0.003242
##
## Model df: 3. Total lags used: 87 Pronóstico
pred <- predict(arimax, newxreg = muestra_exo2)
pred ## $pred
## Qtr1 Qtr2 Qtr3 Qtr4
## 2023 0.12856176 0.07404243 0.02985129 0.01649086
##
## $se
## Qtr1
## 2023 0.06379442