Modelo SARIMA
Roberto Ryder López
20/04/2023
Los modelos ARIMAX
library(readxl)
library(tseries)
library(stargazer)
library(lmtest)
library(car)
library(strucchange)
library(TSstudio)
library(forecast)
library(urca)# Base de datos
data_sarima <- read_excel("data_sarima.xlsx",
sheet = "data2")
head(data_sarima)Se declara la serie de tiempo.
# Tasa de crecimiento logaritmica del PIB
tly <- ts(data_sarima[ , 1],
start = c(1994, 1),
end = c(2022, 4),
freq = 4)
plot(tly)
# Descomposicion de la serie
ts_decompose(tly)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)
adf.test(tly)##
## Augmented Dickey-Fuller Test
##
## data: tly
## Dickey-Fuller = -3.6428, Lag order = 4, p-value = 0.03249
## alternative hypothesis: stationaryadf <- ur.df(tly,
type = "none",
lags = 0)
summary(adf)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.34968 -0.01077 0.04339 0.10282 1.33446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.3213 0.0689 -4.664 8.5e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2317 on 114 degrees of freedom
## Multiple R-squared: 0.1602, Adjusted R-squared: 0.1529
## F-statistic: 21.75 on 1 and 114 DF, p-value: 8.501e-06
##
##
## Value of test-statistic is: -4.6638
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62# Prueba de Phillips-Perron (PP)
pp.test(tly)##
## Phillips-Perron Unit Root Test
##
## data: tly
## Dickey-Fuller Z(alpha) = -44.43, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarypp <- ur.pp(tly)
summary(pp)##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept
##
##
## Call:
## lm(formula = y ~ y.l1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.40573 -0.04623 0.00501 0.06343 1.27821
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.04833 0.02323 2.081 0.0397 *
## y.l1 0.61706 0.07409 8.329 2.16e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2284 on 113 degrees of freedom
## Multiple R-squared: 0.3804, Adjusted R-squared: 0.3749
## F-statistic: 69.37 on 1 and 113 DF, p-value: 2.157e-13
##
##
## Value of test-statistic, type: Z-alpha is: -43.3295
##
## aux. Z statistics
## Z-tau-mu 2.0676# Prueba KPSS
kpss <- ur.kpss(tly)
summary(kpss)##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 4 lags.
##
## Value of test-statistic is: 0.1661
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739La serie es estacionaria.
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.
# Muestra
estacionaria <- cbind(tly)
muestra_end <- window(estacionaria,
start = c(1994, 1),
end = c(2022,4),
freq = 4)3 Correlograma
Se identifican los componentes p (rezagos) y q (medias móviles) con ayuda del correlograma.
# Correlograma
ts_cor(tly)4 Estimación
sarima <- Arima(muestra_end,
order = c(1, 0, 0),
seasonal = list(order = c(0, 0, 1),
period = 4))
summary(sarima)## Series: muestra_end
## ARIMA(1,0,0)(0,0,1)[4] with non-zero mean
##
## Coefficients:
## ar1 sma1 mean
## 0.8792 -0.9998 0.1246
## s.e. 0.0488 0.3007 0.0121
##
## sigma^2 = 0.02783: log likelihood = 37.99
## AIC=-67.98 AICc=-67.62 BIC=-56.97
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01170147 0.1646546 0.07849327 9.319599 52.83753 0.3014549
## ACF1
## Training set -0.07325235coeftest(sarima)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.879188 0.048762 18.0302 < 2.2e-16 ***
## sma1 -0.999766 0.300747 -3.3243 0.0008865 ***
## intercept 0.124578 0.012119 10.2793 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15 Evaluación
5.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(sarima)
adf.test(residuals)##
## Augmented Dickey-Fuller Test
##
## data: residuals
## Dickey-Fuller = -4.9964, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary5.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(sarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(0,0,1)[4] with non-zero mean
## Q* = 2.3453, df = 6, p-value = 0.8854
##
## Model df: 2. Total lags used: 853.56 # Pronóstico
predsarima <- forecast(sarima, h = 8)
summary(predsarima)##
## Forecast method: ARIMA(1,0,0)(0,0,1)[4] with non-zero mean
##
## Model Information:
## Series: muestra_end
## ARIMA(1,0,0)(0,0,1)[4] with non-zero mean
##
## Coefficients:
## ar1 sma1 mean
## 0.8792 -0.9998 0.1246
## s.e. 0.0488 0.3007 0.0121
##
## sigma^2 = 0.02783: log likelihood = 37.99
## AIC=-67.98 AICc=-67.62 BIC=-56.97
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01170147 0.1646546 0.07849327 9.319599 52.83753 0.3014549
## ACF1
## Training set -0.07325235
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2023 Q1 0.2263543 0.009083308 0.4436254 -0.1059330 0.5586417
## 2023 Q2 0.1485613 -0.140582562 0.4377052 -0.2936461 0.5907688
## 2023 Q3 0.1993292 -0.134844329 0.5335027 -0.3117451 0.7104035
## 2023 Q4 0.2436177 -0.121459690 0.6086951 -0.3147200 0.8019553
## 2024 Q1 0.2292362 -0.144676975 0.6031495 -0.3426147 0.8010872
## 2024 Q2 0.2165923 -0.164010278 0.5971948 -0.3654891 0.7986736
## 2024 Q3 0.2054758 -0.180217873 0.5911695 -0.3843918 0.7953434
## 2024 Q4 0.1957024 -0.193881038 0.5852858 -0.4001140 0.7915188