Datos PIB, tomando los primeros 10 años (2014-2023)

Se proceden a cargar las librerías en R

library(readxl)
library(forecast)
library(tseries)
library(urca)
library(nortest)
library(moments)

Leer el archivo de PIB

x= PIB <- read_excel("PIB.xlsx")
Paso 1. Convertir a objeto de serie de Tiempo en R
x=ts(x, start = c(1, 1), end = c(10, 4), frequency = 4)
x
##        Qtr1     Qtr2     Qtr3     Qtr4
## 1  185962.6 190290.0 196779.6 208556.7
## 2  191216.5 196798.2 203824.3 212853.0
## 3  195758.7 200982.0 207001.4 217746.9
## 4  197840.4 203513.3 210525.5 220776.8
## 5  201119.0 209139.1 216587.4 227162.6
## 6  208067.7 215532.9 223374.8 234248.6
## 7  209046.7 179359.4 203148.9 226345.0
## 8  212016.1 212786.4 230209.7 251230.9
## 9  229390.2 238991.9 247215.4 256700.5
## 10 185962.6 190290.0 196779.6 208556.7
class(x)
## [1] "ts"
plot(x, main="Serie de Tiempo", ylab="TRM", col="red")

Gráficas de ACF y PACF de PIB en sus datos originales
acf(ts(x, frequency=1))

pacf(ts(x, frequency = 1))

Paso 2. Prueba DickeyFuller
df1=adf.test(x, k=0) 
df1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  x
## Dickey-Fuller = -3.921, Lag order = 0, p-value = 0.02286
## alternative hypothesis: stationary

Dado que el p-valor es menor que 0.05, se rechaza la hipótesis nula de que la serie temporal tiene una raíz unitaria (no estacionaria) y se acepta la hipótesis alternativa de que la serie es estacionaria.

Paso 3. Modelos

Modelo 1 Arima (2,0,2) —> En sus datos en originales

modelo1=arima(x, order=c(2,0,2))
summary(modelo1)
## 
## Call:
## arima(x = x, order = c(2, 0, 2))
## 
## Coefficients:
##           ar1      ar2     ma1    ma2   intercept
##       -0.1628  -0.5705  0.8663  1.000  211046.807
## s.e.   0.1417   0.1468  0.1365  0.259    3319.213
## 
## sigma^2 estimated as 162741580:  log likelihood = -437.61,  aic = 887.22
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 50.82985 12757.02 9942.656 -0.3575902 4.724756 0.7283658 0.1055697
tsdiag(modelo1)

residuos_1=residuals(modelo1)

Test de LJung-Box

Box.test(residuals(modelo1), type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo1)
## X-squared = 0.48009, df = 1, p-value = 0.4884

La serie tiene ruido blanco, es Estacionaria.

Modelo 2 Arima (2,0,0) = AR(2)

modelo2=arima(x, order=c(2,2,0))
summary(modelo2)
## 
## Call:
## arima(x = x, order = c(2, 2, 0))
## 
## Coefficients:
##           ar1      ar2
##       -0.5723  -0.3701
## s.e.   0.1482   0.1450
## 
## sigma^2 estimated as 5.24e+08:  log likelihood = -435.62,  aic = 877.25
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE     MASE       ACF1
## Training set 185.4796 22310.87 16726.04 -0.2140999 8.111842 1.225294 -0.2675739
tsdiag(modelo2)

Test de LJung-Box

Box.test(residuals(modelo2), type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo2)
## X-squared = 3.0841, df = 1, p-value = 0.07906

La serie tiene ruido blanco, es Estacionaria.

Modelo 3 Arima (1,0,1) = AR(1)

modelo3=arima(x, order=c(1,0,1))
summary(modelo3)
## 
## Call:
## arima(x = x, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1     ma1   intercept
##       0.3294  0.2394  210907.049
## s.e.  0.2375  0.2219    4411.714
## 
## sigma^2 estimated as 234721058:  log likelihood = -442.4,  aic = 892.8
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 297.4338 15320.61 11670.22 -0.3746788 5.528381 0.8549211
##                      ACF1
## Training set -0.001396791
tsdiag(modelo3)

Test de LJung-Box

Box.test(residuals(modelo3), type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo3)
## X-squared = 8.4044e-05, df = 1, p-value = 0.9927

La serie tiene ruido blanco, es Estacionaria.

En resumen: Los 3 modelos pasan, no obstante el mejor modelo es el 1, ya que cuenta con el menor RMSE y AIC.

Pronósticos Arima para el año 11 en trimestre con el modelo1 (2,0,2)
pronostico=forecast::forecast(modelo1, h=4)
pronostico
##       Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 11 Q1       202198.7 185467.3 218930.1 176610.2 227787.1
## 11 Q2       204307.1 184049.3 224565.0 173325.5 235288.8
## 11 Q3       217191.8 196428.4 237955.2 185436.9 248946.7
## 11 Q4       213891.1 191896.1 235886.2 180252.6 247529.7
plot(pronostico)

Modelo Arima (2,0,2) con intercepto
modelo_int=Arima(x, order=c(2,0,2), include.drift = TRUE)
summary(modelo_int)
## Series: x 
## ARIMA(2,0,2) with drift 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2   intercept     drift
##       -0.2174  -0.6478  0.8443  1.0000  198211.688  628.1763
## s.e.   0.1340   0.1388  0.1044  0.1596    5795.436  245.8961
## 
## sigma^2 = 167222815:  log likelihood = -434.82
## AIC=883.65   AICc=887.15   BIC=895.47
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -5.182244 11922.22 9254.743 -0.3399419 4.423942 0.6532451
##                    ACF1
## Training set 0.06667182
tsdiag(modelo_int)

Box.test(residuals(modelo_int), type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo_int)
## X-squared = 0.19148, df = 1, p-value = 0.6617
residuos_int=residuals(modelo_int)
plot(residuos_int)

Se selecciona el modelo_int, ya que cuenta con un menor AIC y RMSE

Pronósticos Arima para el año 11 en trimestre con el modelo_int
pronostico1=forecast::forecast(modelo_int, h=4)
pronostico1
##       Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 11 Q1       210597.3 193642.6 227552.0 184667.3 236527.3
## 11 Q2       220685.5 200849.9 240521.1 190349.5 251021.4
## 11 Q3       234734.1 214683.3 254784.9 204069.1 265399.1
## 11 Q4       226316.8 204975.2 247658.4 193677.6 258956.0
plot(pronostico1)

Mejor Modelo usando la función auto.arima

auto.arima(x, seasonal = F, ic="aic")
## Series: x 
## ARIMA(0,1,0) 
## 
## sigma^2 = 318881474:  log likelihood = -437.15
## AIC=876.31   AICc=876.42   BIC=877.97
modelo5=auto.arima(x, seasonal = F, ic="aic")
summary(modelo5)
## Series: x 
## ARIMA(0,1,0) 
## 
## sigma^2 = 318881474:  log likelihood = -437.15
## AIC=876.31   AICc=876.42   BIC=877.97
## 
## Training set error measures:
##                    ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 569.5014 17632.62 13314.02 -0.06582525 6.421169 0.9397687
##                    ACF1
## Training set -0.1539005
residuos_aut=resid(modelo5)

En este caso el mejor modelo sigue siendo el modelo_int, porque tienen el menor RMSE; no obstante, para el criterio AIC el mejor modelo es el de la función auto.arima ya que cuenta con el menor AIC.

Pronósticos Arima para el año 11 en trimestres con el modelo auto.arima criterio AIC
pronostico2=forecast::forecast(modelo5, h=4)
pronostico2
##       Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 11 Q1       208556.7 185671.7 231441.7 173557.2 243556.3
## 11 Q2       208556.7 176192.5 240921.0 159059.9 258053.6
## 11 Q3       208556.7 168918.8 248194.7 147935.7 269177.8
## 11 Q4       208556.7 162786.8 254326.7 138557.6 278555.9
plot(pronostico2)

Paso 4. Determinando si la serie es Estacional

modelo6=auto.arima(x, seasonal = T, ic="aic")
summary(modelo6)
## Series: x 
## ARIMA(2,0,1)(0,1,0)[4] 
## 
## Coefficients:
##          ar1      ar2      ma1
##       1.6859  -0.8802  -0.8094
## s.e.  0.0944   0.0814   0.1283
## 
## sigma^2 = 124061368:  log likelihood = -386.26
## AIC=780.51   AICc=781.8   BIC=786.84
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 1638.512 10116.85 6137.386 0.5990931 2.974606 0.4332068 -0.0724326

La serie es Estacionaria, pero no es Estacional, la parte de medias móviles están en 0, adicional el modelo 6 cuenta con un menor AIC y RMSE

Pronósticos Arima para el año 11 en trimestres con el modelo 6
pronostico3=forecast(modelo6, h=4, level=95)
pronostico3
##       Point Forecast    Lo 95    Hi 95
## 11 Q1       157674.4 135843.8 179505.0
## 11 Q2       184974.3 155944.4 214004.2
## 11 Q3       212717.3 180890.9 244543.7
## 11 Q4       240105.9 207865.4 272346.4
autoplot(pronostico3)

Prueba de Normalidad para los mejores modelos

Modelo 1: Tiene el menor AIC y RMSE

qqnorm(residuos_1)
qqline(residuos_1)

hist(residuos_1)

plot((density(residuos_1)), col="blue", las=1, lwd=4)

jarque.bera.test(residuos_1)
## 
##  Jarque Bera Test
## 
## data:  residuos_1
## X-squared = 1.739, df = 2, p-value = 0.4192
shapiro.test(residuos_1)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_1
## W = 0.97825, p-value = 0.6247
lillie.test(residuos_1)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos_1
## D = 0.076139, p-value = 0.8095

Análisis:

  1. se observa que los puntos del gráfico se alinean aproximadamente en una línea recta, esto indica que los datos siguen la distribución teórica, ya que cuanto más cerca estén los puntos de la línea negra, mejor será el ajuste; no obstante,hay algunos puntos que se desvían sistemáticamente de la línea

  2. El la gráfica del histograma se observa que hay un ligero sesgo hacia la derecha, sin embargo, los datos no son totalmente asimétricos esto se observa más detalladamente con la gráfica de densidad, donde se logra observar una forma de campana.

  3. Jarque-Bera Test: Dado que el p-valor es mayor que el nivel de significancia se acepta la hipótesis nula de que los residuos tienen una distribución normal.

  4. Shapiro-Wilk Normality Test: Dado que el estadístico de la prueba de Shapiro-Wilk (W) es un valor cercano a 1 indica que los datos son aproximadamente normales, adicional el p-value es mayor al nivel de significancia se acepta la hipótesis nula de que los residuos son normalmente distribuidos.

  5. Lilliefors (Kolmogorov-Smirnov) Normality Test: Se acpeta la hipótesis nula de que los residuos son normalmente distribuidos,ya que el P-value el mayor al nivel de significancia.

Moldelo Arima con intercepto —> tiene menor AIC y RMSE

qqnorm(residuos_int)
qqline(residuos_int)

hist(residuos_int)

plot((density(residuos_int)), col="blue", las=1, lwd=4)

jarque.bera.test(residuos_int)
## 
##  Jarque Bera Test
## 
## data:  residuos_int
## X-squared = 1.4804, df = 2, p-value = 0.477
shapiro.test(residuos_int)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_int
## W = 0.9554, p-value = 0.1164
lillie.test(residuos_int)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos_int
## D = 0.13953, p-value = 0.0482

Análisis:

  1. se observa que los puntos del gráfico se alinean aproximadamente en una línea recta; no obstante, hay algunos puntos que se desvían sistemáticamente de la línea

  2. El la gráfica del histograma se observa que hay un ligero sesgo hacia la izquierda, sin embargo los datos no son totalmente asimétricos esto se observa más detalladamente con la gráfica de densidad, donde se logra observar una forma de campana.

  3. Jarque-Bera Test: Se acepta la hipótesis nula de que los residuos tienen una distribución normal.

  4. Shapiro-Wilk Normality Test: Dado que el estadístico de la prueba de Shapiro-Wilk (W) es un valor cercano a 1y el p-value es mayor al nivel de significancia se acepta la hipótesis nula de que los residuos son normalmente distribuidos.

  5. Lilliefors (Kolmogorov-Smirnov) Normality Test: Dado que el p-valor es menor que el nivel de significancia, se puede rechazar la hipótesis nula. Esto sugiere que hay evidencia suficiente para concluir que los residuos no son normalmente distribuidos.

Moldelo 5 auto.arima —> tiene menor AIC

qqnorm(residuos_aut)
qqline(residuos_aut)

hist(residuos_aut)

plot((density(residuos_aut)), col="blue", las=1, lwd=4)

jarque.bera.test(residuos_aut)
## 
##  Jarque Bera Test
## 
## data:  residuos_aut
## X-squared = 57.314, df = 2, p-value = 3.584e-13
shapiro.test(residuos_aut)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_aut
## W = 0.79199, p-value = 4.632e-06
lillie.test(residuos_aut)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos_aut
## D = 0.28338, p-value = 9.571e-09

Análisis:

Con el modelo 5 auto.arima no se cumplen el supuesto de normalidad, se observa en las gráficas la presencia de valores atípicos, los cuales se encuentran muy alejados del resto de datos, es decir, no sigue una distribución normal. El test Jarque-Bera, Shapiro-Wilk y Lilliefors, rechazan la hipótesis nula, hay evidencia sufiente para aceptar la hipótesis alternativa de que los residuos no son normalmente distribuidos.

Conclusiones

Estacionariedad de la serie (PIB): Según la prueba Dickey-Fuller la serie es Estacionaria, dado que el p-valor = 0.02286 es menor al 0.05. Esto implica que las características estadísticas de la serie no varían con el tiempo, lo que facilita el modelado de la serie sin necesidad de diferenciación.

Selección del mejor modelo: Entre los diferentes modelos ARIMA que se ejecutaron, el modelo ARIMA(2,0,2) con intercepto es el mejor, ya que cuenta con un menor RMSE (11922.22) y AIC (883.65). Estos indicadores sugieren que este modelo ofrece un ajuste más preciso y eficiente en comparación con los demás.

Estacionalidad: Aunque la serie es estacionaria, el modelo estacional ARIMA(2,0,1)(0,1,0)[4] muestra que la serie no es estacional, ya que no presenta patrones repetitivos a intervalos regulares. Las medias móviles estimadas están cercanas a 0, lo que refuerza la conclusión de que no hay un componente estacional significativo en la serie.

Resultados de las pruebas de normalidad: El modelo ARIMA(2,0,2) con intercepto muestra un comportamiento aceptable en términos de normalidad en los residuos, con resultados favorables en las pruebas de Jarque-Bera y Shapiro-Wilk; no obstante,el modelo obtenido con la función auto.arima no pasa las pruebas de normalidad, con p-valores muy bajos en las pruebas de Jarque-Bera, Shapiro-Wilk y Lilliefors, lo que indica que los residuos no siguen una distribución normal.

Por último, el modelo obtenido con auto.arima tiene el menor AIC; no obstante, el modelo ARIMA(2,0,2) con intercepto es más adecuado para la serie temporal del PIB debido a su menor RMSE, lo que indica un mejor ajuste de los datos. Además, sus residuos muestran un comportamiento más cercano a la normalidad, tal y como se puedo apreciar en la gráficas y test