Serie de Tiempo Económica

Para el desarrollo de la actividad utilizaremos la serie económica “Reservas Internacionales Brutas - Banco de La República” la cual fue utilizada en la primera y segunda actividad.

Serie de tiempo

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug     Sep
## 2014                                                 46099.3 46712.1 47115.1
## 2015 47072.6 47060.3 46920.0 47179.9 47012.7 46981.6 46852.8 46710.2 46732.6
## 2016 46740.5 46852.9 47228.6 47297.8 47538.5 47029.6 47096.7 47043.1 47113.0
## 2017 46989.3 46992.3 46937.4 46941.8 47159.3 47242.5 47540.7 47577.2 47524.6
## 2018 47791.1 47642.5 47614.4 47510.7 47607.0 47497.2 47526.2 47543.5 47520.0
## 2019 49217.2 50503.4 51267.3 51528.3 51978.8 52449.3 52372.4 52998.6 52874.6
## 2020 53469.7 53681.6 53340.7 53868.1 56367.6 56629.2 56978.8 57193.1 56987.1
## 2021 59010.3 58982.0 58908.7 59094.6 59152.9 58925.5 58885.6 58854.4 58730.2
## 2022 58272.9 58289.2 58010.3 57480.5 57634.7 57170.9 57359.7 56989.5 56339.9
## 2023 57802.6 57388.6 57990.1 58058.0 57721.3 57866.2 58152.5 57961.8 57595.4
## 2024 59697.1 59480.3 60023.6 59847.0 60600.5 60931.5                        
##          Oct     Nov     Dec
## 2014 47373.3 47388.6 47328.1
## 2015 46833.8 46766.4 46740.4
## 2016 46974.3 46751.1 46682.8
## 2017 47416.9 47403.9 47637.2
## 2018 47505.6 47761.2 48401.5
## 2019 53090.5 52957.9 53174.2
## 2020 56928.1 57254.4 59039.3
## 2021 58800.2 58540.2 58587.8
## 2022 56398.9 56991.7 57290.1
## 2023 57488.2 58609.4 59639.2
## 2024

Estadísticos Descriptivos

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   46099   47239   52411   52468   57742   60932

Gráfico de la Serie de Tiempo

# Convertir la serie de tiempo a un dataframe
df_ts <- data.frame(
  fecha = seq(from = as.Date("2014-07-01"), by = "month", length.out = length(reserv_intern_brut_col)),
  valor = as.numeric(reserv_intern_brut_col)
)

ggplot(data = df_ts, aes(x = fecha, y = valor)) +
  geom_line(color = "blue") +  # Línea para la serie de tiempo
  geom_point(color = "blue") +   # Puntos en cada observación
  labs(title = "Reservas Internacionales Brutas Banco de la Republica",
       x = "Fecha",
       y = "Valor") +
  theme_minimal()  # Tema sencillo

Se observa una tendencia creciente en la gráfica de la serie de tiempo lo cual es un indicio de que la serie no es estacionaria.

Funciones de Autocorrelación y Autocorrelación Parcial - Datos Originales

acf(ts(reserv_intern_brut_col, frequency = 1), main = "Función de autocorrelación - Reservas Internacionales")

pacf(ts(reserv_intern_brut_col, frequency = 1), main = "Función de autocorrelación Parcial - Reservas Internacionales")

ndiffs(reserv_intern_brut_col)
## [1] 1

Prueba de Dickey-Fuller Datos Originales

dick_full_1 <- adf.test(reserv_intern_brut_col, alternative = "stationary", k=0)
dick_full_1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  reserv_intern_brut_col
## Dickey-Fuller = -1.3819, Lag order = 0, p-value = 0.8333
## alternative hypothesis: stationary

La serie tiene raiz unitaria, es decir no es estacionaria. el valor p es mayor a 0.05, por lo que no se puede rechazar la hipótesis nula y, por lo tanto, la serie es no estacionaria.

Diferencias de la Serie

serie_diff <- diff(reserv_intern_brut_col)
serie_diff
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2014                                                   612.8  403.0  258.2
## 2015 -255.5  -12.3 -140.3  259.9 -167.2  -31.1 -128.8 -142.6   22.4  101.2
## 2016    0.1  112.4  375.7   69.2  240.7 -508.9   67.1  -53.6   69.9 -138.7
## 2017  306.5    3.0  -54.9    4.4  217.5   83.2  298.2   36.5  -52.6 -107.7
## 2018  153.9 -148.6  -28.1 -103.7   96.3 -109.8   29.0   17.3  -23.5  -14.4
## 2019  815.7 1286.2  763.9  261.0  450.5  470.5  -76.9  626.2 -124.0  215.9
## 2020  295.5  211.9 -340.9  527.4 2499.5  261.6  349.6  214.3 -206.0  -59.0
## 2021  -29.0  -28.3  -73.3  185.9   58.3 -227.4  -39.9  -31.2 -124.2   70.0
## 2022 -314.9   16.3 -278.9 -529.8  154.2 -463.8  188.8 -370.2 -649.6   59.0
## 2023  512.5 -414.0  601.5   67.9 -336.7  144.9  286.3 -190.7 -366.4 -107.2
## 2024   57.9 -216.8  543.3 -176.6  753.5  331.0                            
##         Nov    Dec
## 2014   15.3  -60.5
## 2015  -67.4  -26.0
## 2016 -223.2  -68.3
## 2017  -13.0  233.3
## 2018  255.6  640.3
## 2019 -132.6  216.3
## 2020  326.3 1784.9
## 2021 -260.0   47.6
## 2022  592.8  298.4
## 2023 1121.2 1029.8
## 2024

La diferencia de la serie se hace a 1 periodo

Gráfico de las diferencias de la Serie de Tiempo

df_ts_diff <- data.frame(
  fecha = seq(from = as.Date("2014-07-08"), by = "month", length.out = length(serie_diff)),
  valor = as.numeric(serie_diff)
)

ggplot(data = df_ts_diff, aes(x = fecha, y = valor)) +
  geom_line(color = "blue") +  # Línea para la serie de tiempo
  geom_point(color = "blue") +   # Puntos en cada observación
  labs(title = "Reservas Internacionales Brutas Banco de la Republica - Diferencias",
       x = "Fecha",
       y = "Valor") +
  theme_minimal()  # Tema sencillo

Funciones de Autocorrelación y Autocorrelación Parcial - Datos Diferenciados

acf(ts(serie_diff, frequency = 1))

pacf(ts(serie_diff, frequency = 1))

ndiffs(serie_diff)
## [1] 0

#La serie ya es estacionaria porque se necesitan 0 diferencias.

Prueba de Dickey-Fuller Datos Diferenciados

dick_full_2 <- adf.test(serie_diff, alternative = "stationary", k=0)
## Warning in adf.test(serie_diff, alternative = "stationary", k = 0): p-value
## smaller than printed p-value
dick_full_2
## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie_diff
## Dickey-Fuller = -8.5411, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

Modelos Arima

Modelo 1

modelo1 <- arima(reserv_intern_brut_col, order=c(1, 1, 2))
summary(modelo1)
## 
## Call:
## arima(x = reserv_intern_brut_col, order = c(1, 1, 2))
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.9457  -0.6980  -0.1453
## s.e.  0.0623   0.1129   0.0962
## 
## sigma^2 estimated as 170580:  log likelihood = -885.78,  aic = 1779.56
## 
## Training set error measures:
##                   ME     RMSE    MAE        MPE      MAPE      MASE        ACF1
## Training set 48.6745 411.3106 273.12 0.09149361 0.5056718 0.9991108 -0.01609318
tsdiag(modelo1)

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

Modelo 2

modelo2 <- Arima(reserv_intern_brut_col, order=c(1, 1, 1), include.drift = TRUE)
summary(modelo2)
## Series: reserv_intern_brut_col 
## ARIMA(1,1,1) with drift 
## 
## Coefficients:
##          ar1      ma1     drift
##       0.2718  -0.0376  126.6106
## s.e.  0.4219   0.4398   49.5737
## 
## sigma^2 = 172425:  log likelihood = -884.79
## AIC=1777.58   AICc=1777.93   BIC=1788.7
## 
## Training set error measures:
##                      ME     RMSE      MAE          MPE      MAPE     MASE
## Training set -0.7578052 408.2611 276.0999 -0.007492145 0.5136575 0.154017
##                      ACF1
## Training set 0.0005664837
tsdiag(modelo2)

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

Modelo 3

modelo3 <- arima(reserv_intern_brut_col, order=c(0, 0, 1))
summary(modelo3)
## 
## Call:
## arima(x = reserv_intern_brut_col, order = c(0, 0, 1))
## 
## Coefficients:
##          ma1   intercept
##       1.0000  52467.6687
## s.e.  0.0218    481.0808
## 
## sigma^2 estimated as 7000526:  log likelihood = -1118.36,  aic = 2242.72
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE     MASE     ACF1
## Training set 17.08639 2645.851 2457.626 -0.4570596 4.710502 8.990336 0.919613
tsdiag(modelo3)

Box.test(residuals(modelo3), type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo3)
## X-squared = 104.04, df = 1, p-value < 2.2e-16

Modelo 4

modelo4 <- Arima(reserv_intern_brut_col, order=c(0, 1, 2), include.drift = TRUE)
summary(modelo4)
## Series: reserv_intern_brut_col 
## ARIMA(0,1,2) with drift 
## 
## Coefficients:
##          ma1     ma2     drift
##       0.2358  0.0589  126.3357
## s.e.  0.0920  0.0904   48.5615
## 
## sigma^2 = 172449:  log likelihood = -884.8
## AIC=1777.6   AICc=1777.95   BIC=1788.72
## 
## Training set error measures:
##                      ME     RMSE      MAE          MPE      MAPE      MASE
## Training set -0.6338815 408.2902 276.3929 -0.007391297 0.5142357 0.1541805
##                      ACF1
## Training set 0.0001024154
tsdiag(modelo4)

Box.test(residuals(modelo4), type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo4)
## X-squared = 1.2904e-06, df = 1, p-value = 0.9991

Modelo 5

modelo5 <- auto.arima(reserv_intern_brut_col, seasonal = FALSE, ic = "aic")
summary(modelo5)
## Series: reserv_intern_brut_col 
## ARIMA(1,1,0) with drift 
## 
## Coefficients:
##          ar1     drift
##       0.2362  126.4354
## s.e.  0.0892   49.0828
## 
## sigma^2 = 170963:  log likelihood = -884.8
## AIC=1775.59   AICc=1775.8   BIC=1783.93
## 
## Training set error measures:
##                      ME     RMSE      MAE          MPE      MAPE      MASE
## Training set -0.6887114 408.2759 276.3436 -0.007415118 0.5141781 0.1541529
##                      ACF1
## Training set -0.001267992
tsdiag(modelo5)

Box.test(residuals(modelo5), type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo5)
## X-squared = 0.0001978, df = 1, p-value = 0.9888

Modelo 6 - Estacional

modelo6 <- auto.arima(reserv_intern_brut_col, seasonal = TRUE, ic = "aic")
summary(modelo6)
## Series: reserv_intern_brut_col 
## ARIMA(1,1,0) with drift 
## 
## Coefficients:
##          ar1     drift
##       0.2362  126.4354
## s.e.  0.0892   49.0828
## 
## sigma^2 = 170963:  log likelihood = -884.8
## AIC=1775.59   AICc=1775.8   BIC=1783.93
## 
## Training set error measures:
##                      ME     RMSE      MAE          MPE      MAPE      MASE
## Training set -0.6887114 408.2759 276.3436 -0.007415118 0.5141781 0.1541529
##                      ACF1
## Training set -0.001267992
tsdiag(modelo6)

Box.test(residuals(modelo6), type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo6)
## X-squared = 0.0001978, df = 1, p-value = 0.9888

Resumen de los modelos AIC, RMSE y BIC

## AIC Modelo 1: 1779.563
## AIC Modelo 2: 1777.582
## AIC Modelo 3: 2242.721
## AIC Modelo 4: 1777.599
## AIC Modelo 5: 1775.59  Mejor Modelo Criterio AIC
## BIC Modelo 1: 1790.68
## BIC Modelo 2: 1788.699
## BIC Modelo 3: 2251.083
## BIC Modelo 4: 1788.716
## BIC Modelo 5: 1783.928
## RMSE Modelo 1: 411.3106
## RMSE Modelo 2: 408.2611 Mejor Modelo Criterio RMSE
## RMSE Modelo 3: 410.7178
## RMSE Modelo 4: 414.9445
## RMSE Modelo 5: 408.2759

Errores modelo seleccionado

erroresmod5 <- residuals(modelo5)
plot(erroresmod5, main = "Errores del Modelo ARIMA No. 5 - Mejor AIC")

erroresmod2 <- residuals(modelo2) 
plot(erroresmod2, main = "Errores del Modelo ARIMA No. 2 - Mejor RMSE")

Pruebas de normalidad del modelo con mejor AIC (Modelo 5)

Gráfico Q-Q

qqnorm(erroresmod5, main = "Gráfico Q-Q de los Residuos - Modelo ARIMA 5")
qqline(erroresmod5, col = "red")

Coeficiente de asimetría

coef_asimetria_5 <- skewness(erroresmod5)
cat("Coeficiente de Asimetría Modelo ARIMA 5:", coef_asimetria_5, "\n")
## Coeficiente de Asimetría Modelo ARIMA 5: 2.125571

Coeficiente de curtosis

coef_curtosis_5 <- kurtosis(erroresmod5)
cat("Coeficiente de Curtosis Modelo ARIMA 5:", coef_curtosis_5, "\n")
## Coeficiente de Curtosis Modelo ARIMA 5: 8.324607

Prueba de Shapiro-Wilk

shapiro_test_5 <- shapiro.test(erroresmod5)
print(shapiro_test_5)
## 
##  Shapiro-Wilk normality test
## 
## data:  erroresmod5
## W = 0.84605, p-value = 7.742e-10

Prueba de Kolmogorov-Smirnov-Lilliefors

lilliefors_test_5 <- lillie.test(erroresmod5)
print(lilliefors_test_5)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  erroresmod5
## D = 0.13572, p-value = 1.149e-05

Pruebas de normalidad del modelo con mejor RMSE (Modelo 2)

Gráfico Q-Q

qqnorm(erroresmod2, main = "Gráfico Q-Q de los Residuos - Modelo ARIMA 2")
qqline(erroresmod2, col = "red")

Coeficiente de asimetría

coef_asimetria_2 <- skewness(erroresmod2)
cat("Coeficiente de Asimetría Modelo ARIMA 2:", coef_asimetria_2, "\n")
## Coeficiente de Asimetría Modelo ARIMA 2: 2.137306

Coeficiente de curtosis

coef_curtosis_2 <- kurtosis(erroresmod2)
cat("Coeficiente de Curtosis Modelo ARIMA 2:", coef_curtosis_2, "\n")
## Coeficiente de Curtosis Modelo ARIMA 2: 8.401693

Prueba de Shapiro-Wilk

shapiro_test_2 <- shapiro.test(erroresmod2)
print(shapiro_test_2)
## 
##  Shapiro-Wilk normality test
## 
## data:  erroresmod2
## W = 0.84482, p-value = 6.931e-10

Prueba de Kolmogorov-Smirnov-Lilliefors

lilliefors_test_2 <- lillie.test(erroresmod2)
print(lilliefors_test_2)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  erroresmod2
## D = 0.13681, p-value = 9.233e-06

Pronostico año siguiente utilizando modelo con mejor AIC

# Predecir los siguientes 12 períodos
pronostico_proximo_año <- forecast(modelo5, h = 12)
pronostico_proximo_año
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jul 2024       61106.25 60576.36 61636.14 60295.85 61916.65
## Aug 2024       61244.10 60401.56 62086.63 59955.55 62532.64
## Sep 2024       61373.23 60287.62 62458.83 59712.93 63033.52
## Oct 2024       61500.30 60213.12 62787.48 59531.73 63468.87
## Nov 2024       61626.88 60164.90 63088.87 59390.97 63862.80
## Dec 2024       61753.36 60135.17 63371.54 59278.56 64228.15
## Jan 2025       61879.80 60119.19 63640.41 59187.18 64572.42
## Feb 2025       62006.24 60113.88 63898.59 59112.13 64900.35
## Mar 2025       62132.67 60117.16 64148.18 59050.21 65215.13
## Apr 2025       62259.11 60127.54 64390.67 58999.16 65519.05
## May 2025       62385.54 60143.93 64627.16 58957.29 65813.80
## Jun 2025       62511.98 60165.47 64858.49 58923.30 66100.66
plot(pronostico_proximo_año)

Evaluar rendimiento y precision de las predicciones

# Ajuste del modelo
accuracy(modelo4)
##                      ME     RMSE      MAE          MPE      MAPE      MASE
## Training set -0.6338815 408.2902 276.3929 -0.007391297 0.5142357 0.1541805
##                      ACF1
## Training set 0.0001024154

CONCLUSIONES

  • Todos los modelos tienen ruido blanco conforme al Box-Ljung test aplicado a cada uno de estos.
  • La serie no es estacional considerando los modelos auto-arima con estacionalidad y no estacionalidad. Los modelos generan el mismo RMSE y el mismo AIC.
  • El buen ajuste del modelo seleccionado permite hacer predicciones con intervalos de confianza no tan amplios lo que permite reducir la incertidumbre en la predicción.
  • El coeficiente de asimetria del modelo seleccionado nos indica que los residuos del modelo ARIMA tienen una distribución sesgada hacia la derecha, lo que significa que hay una cola más larga en el lado derecho de la distribución de los residuos. Este valor positivo sugiere que hay algunas observaciones de residuos más extremas en el lado derecho.
  • En cuanto a la prueba Shapiro-Wilk el valor de W es 0.84605, que es menor que 1. Esto sugiere que hay una desviación importante de la normalidad en los datos.
  • El valor p de la prueba shapiro wilk es muy pequeña lo cual indica que es extremadamente improbable que los datos provengan de una distribución normal si la hipótesis nula fuera verdadera.