# Cargar las bibliotecas necesarias
library(tseries)## Warning: package 'tseries' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3library(forecast)## Warning: package 'forecast' was built under R version 4.3.31. Captura los datos
t <- 1:90
yt <- c(235.000, 239.000, 244.090, 252.731, 264.377, 277.934, 286.687, 295.629, 310.444, 325.112,
336.291, 344.459, 355.399, 367.691, 384.003, 398.042, 412.969, 422.901, 434.960, 445.853,
455.853, 465.584, 477.894, 491.408, 507.712, 517.237, 524.349, 532.104, 538.097, 544.948,
551.925, 557.929, 564.285, 572.164, 582.926, 595.295, 607.028, 617.541, 622.941, 633.436,
647.371, 658.230, 670.777, 685.457, 690.992, 693.557, 700.675, 712.710, 726.513, 736.429,
743.203, 751.227, 764.265, 777.852, 791.070, 805.844, 815.122, 822.905, 830.663, 839.600,
846.962, 853.830, 860.840, 871.075, 877.792, 881.143, 884.226, 890.208, 894.966, 901.288,
913.138, 922.511, 930.786, 941.306, 950.305, 952.373, 960.042, 968.100, 972.477, 977.408,
977.602, 979.505, 982.934, 985.833, 991.350, 996.291, 1003.100, 1010.320, 1018.420, 1029.480)
data <- data.frame(t, yt)
data## t yt
## 1 1 235.000
## 2 2 239.000
## 3 3 244.090
## 4 4 252.731
## 5 5 264.377
## 6 6 277.934
## 7 7 286.687
## 8 8 295.629
## 9 9 310.444
## 10 10 325.112
## 11 11 336.291
## 12 12 344.459
## 13 13 355.399
## 14 14 367.691
## 15 15 384.003
## 16 16 398.042
## 17 17 412.969
## 18 18 422.901
## 19 19 434.960
## 20 20 445.853
## 21 21 455.853
## 22 22 465.584
## 23 23 477.894
## 24 24 491.408
## 25 25 507.712
## 26 26 517.237
## 27 27 524.349
## 28 28 532.104
## 29 29 538.097
## 30 30 544.948
## 31 31 551.925
## 32 32 557.929
## 33 33 564.285
## 34 34 572.164
## 35 35 582.926
## 36 36 595.295
## 37 37 607.028
## 38 38 617.541
## 39 39 622.941
## 40 40 633.436
## 41 41 647.371
## 42 42 658.230
## 43 43 670.777
## 44 44 685.457
## 45 45 690.992
## 46 46 693.557
## 47 47 700.675
## 48 48 712.710
## 49 49 726.513
## 50 50 736.429
## 51 51 743.203
## 52 52 751.227
## 53 53 764.265
## 54 54 777.852
## 55 55 791.070
## 56 56 805.844
## 57 57 815.122
## 58 58 822.905
## 59 59 830.663
## 60 60 839.600
## 61 61 846.962
## 62 62 853.830
## 63 63 860.840
## 64 64 871.075
## 65 65 877.792
## 66 66 881.143
## 67 67 884.226
## 68 68 890.208
## 69 69 894.966
## 70 70 901.288
## 71 71 913.138
## 72 72 922.511
## 73 73 930.786
## 74 74 941.306
## 75 75 950.305
## 76 76 952.373
## 77 77 960.042
## 78 78 968.100
## 79 79 972.477
## 80 80 977.408
## 81 81 977.602
## 82 82 979.505
## 83 83 982.934
## 84 84 985.833
## 85 85 991.350
## 86 86 996.291
## 87 87 1003.100
## 88 88 1010.320
## 89 89 1018.420
## 90 90 1029.4802. Gráfica de los datos Y verificamos si son estacionarios
ggplot(data, aes(x=t, y=yt)) +
geom_line() +
geom_point() +
labs(title="Datos de yt en función de t", x="t", y="yt") +
theme_minimal()
# Prueba de Dickey-Fuller aumentada
adf_test <- adf.test(data$yt)
# Mostramos los resultados de la prueba ADF
print(adf_test)##
## Augmented Dickey-Fuller Test
##
## data: data$yt
## Dickey-Fuller = -0.74411, Lag order = 4, p-value = 0.9632
## alternative hypothesis: stationaryUn p-valor de 0.9632 es mucho mayor que 0.05, lo que indica que no podemos rechazar la hipótesis nula de que los datos tienen una raíz unitaria. Esto sugiere que los datos no son estacionarios.
3. En caso de que no sean estacionarios, realizar una diferencia, realizar una diferencia
yt_diff1 <- diff(data$yt, differences = 1)
# Agregamos los datos diferenciados al dataframe
data_diff1 <- data.frame(t = t[-1], yt_diff1 = yt_diff1)
# Gráficamos de los datos diferenciados (primer orden)
ggplot(data_diff1, aes(x=t, y=yt_diff1)) +
geom_line() +
geom_point() +
labs(title="Datos Diferenciados (Primer Orden) de yt en función de t", x="t", y="yt_diff1") +
theme_minimal()
# Prueba de Dickey-Fuller aumentada en los datos diferenciados (primer orden)
adf_test_diff1 <- adf.test(data_diff1$yt_diff1)
# Mostramos resultados de la prueba ADF en los datos diferenciados (primer orden)
print(adf_test_diff1)##
## Augmented Dickey-Fuller Test
##
## data: data_diff1$yt_diff1
## Dickey-Fuller = -3.3802, Lag order = 4, p-value = 0.06342
## alternative hypothesis: stationaryp-valor: Un p-valor de 0.06342 es ligeramente superior a 0.05, lo que indica que no podemos rechazar la hipótesis nula de que los datos tienen una raíz unitaria al nivel del 5%. Sin embargo, está muy cerca del umbral, lo que podría sugerir que los datos están casi estacionarios o podrían ser considerados estacionarios al nivel del 10%.
realizamos una segunda diferenciaa
yt_diff2 <- diff(data$yt, differences = 2)
# Agregamos los datos diferenciados de segundo orden al dataframe
data_diff2 <- data.frame(t = t[-c(1,2)], yt_diff2 = yt_diff2)
# Gráficamos de los datos diferenciados (segundo orden)
ggplot(data_diff2, aes(x=t, y=yt_diff2)) +
geom_line() +
geom_point() +
labs(title="Datos Diferenciados (Segundo Orden) de yt en función de t", x="t", y="yt_diff2") +
theme_minimal()
# Prueba de Dickey-Fuller aumentada en los datos diferenciados (segundo orden)
adf_test_diff2 <- adf.test(data_diff2$yt_diff2)## Warning in adf.test(data_diff2$yt_diff2): p-value smaller than printed p-value# Mostramos resultados de la prueba ADF en los datos diferenciados (segundo orden)
print(adf_test_diff2)##
## Augmented Dickey-Fuller Test
##
## data: data_diff2$yt_diff2
## Dickey-Fuller = -6.052, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryInterpretación final Primera diferenciación: Los datos están cerca de ser estacionarios. Segunda diferenciación: Los datos son estacionarios, como lo indica el extremadamente bajo p-valor.
4. Identificar un modelo tentativo mediante las funciones de autocorrelación muestral y autocorrelación parcial muestral
# Calculamos y graficamos la función de autocorrelación muestral (ACF)
acf_yt_diff2 <- acf(data_diff2$yt_diff2, main="ACF de la segunda diferencia", ylim=c(-1,1))
pacf_yt_diff2 <- pacf(data_diff2$yt_diff2, main="PACF de la segunda diferencia", ylim=c(-1,1))
5. realizar el ajuste del modelo identificado
# Ajustamos el modelo ARIMA(1,2,2)
modelo_arima <- arima(data$yt, order = c(1,2,2))
# Resumen del modelo para evaluar los parámetros
summary(modelo_arima)##
## Call:
## arima(x = data$yt, order = c(1, 2, 2))
##
## Coefficients:
## ar1 ma1 ma2
## 0.2920 -0.4933 -0.4042
## s.e. 0.1828 0.1658 0.1316
##
## sigma^2 estimated as 7.25: log likelihood = -212.63, aic = 433.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07518209 2.662657 2.172112 0.03860375 0.3796592 0.2433264
## ACF1
## Training set 0.0060996326. Realizar la valoración estadística del modelo, realizando la prueba de hipótesis de los parámetros del modelo, determinar si el modelo es adecuado mediante la prueba de Ljung-Box y graficar los residuales para verificar que el modelo es adecuado a los datos.
# Obtenemos los residuos del modelo ajustado
residuos <- residuals(modelo_arima)
residuos## Time Series:
## Start = 1
## End = 90
## Frequency = 1
## [1] 0.1050952 -0.3063403 0.9652608 3.1762636 3.2901429 3.3058523
## [7] -2.6105078 1.6583446 5.4012777 1.2258653 -0.7389443 -1.8255784
## [13] 2.4635584 0.9775224 5.0302209 -0.6321325 3.2401343 -3.9166013
## [19] 2.9703846 -1.9077810 -0.2869612 -0.9137542 2.0895792 1.1036799
## [25] 3.8148467 -5.2679961 -1.4792345 -1.5035313 -3.2826668 -0.8484330
## [31] -1.8661792 -2.2695407 -1.2344533 -0.1045204 1.8872905 1.6522280
## [37] 0.4713960 -0.1342851 -4.6313912 4.2497539 2.1753196 -1.2901760
## [43] 2.8288417 2.5132779 -7.3843840 -2.9252208 0.9927819 2.8946289
## [49] 2.1610644 -2.1672579 -2.2023300 0.2051475 3.8598579 1.0716711
## [55] 1.5595032 2.8661005 -3.9061845 -0.6583826 -1.4922056 0.1841212
## [61] -2.4316110 -1.1590957 -1.2684235 2.0893147 -3.9418013 -3.4385474
## [67] -2.5746866 0.3172645 -2.9547502 0.5921597 4.1689995 -1.7953019
## [73] 0.4249479 2.0494891 -0.9937836 -6.1486051 4.1901243 -1.6650536
## [79] -2.9221028 -0.4856517 -6.3195371 -0.2214101 -1.6368110 -1.8724675
## [85] 1.1874610 -1.5116060 1.7705689 0.1278890 1.5388052 3.5137952# Realizamos la prueba de Ljung-Box para los residuos
ljung_box_test <- Box.test(residuos, type = "Ljung-Box")
ljung_box_test##
## Box-Ljung test
##
## data: residuos
## X-squared = 0.0034614, df = 1, p-value = 0.9531p-value = 0.9531: Esto es mucho mayor que 0.05, lo que indica que no hay suficiente evidencia para rechazar la hipótesis nula de que no hay autocorrelación en los residuos del modelo. En otras palabras, los residuos del modelo ARIMA(1,2,2) parecen no tener autocorrelación significativa, lo que sugiere que el modelo es adecuado.
# Graficamos los residuos del modelo
ggplot(data.frame(t = 1:length(residuos), residuos = residuos), aes(x = t, y = residuos)) +
geom_line() +
geom_point() +
labs(title = "Residuos del modelo ARIMA(1,2,2)", x = "Tiempo", y = "Residuo") +
theme_minimal()## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# Graficamos la ACF de los residuos
acf(residuos, main = "ACF de los residuos")
# Graficamos la PACF de los residuos
pacf(residuos, main = "PACF de los residuos")
7.Realizar los pronósticos puntuales y los intervalos de predicción al 95%, de los periodos t=91 al t=100.
forecast_arima <- forecast(modelo_arima, h = 10)
print(forecast_arima)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 91 1039.049 1035.598 1042.500 1033.772 1044.326
## 92 1046.762 1039.661 1053.864 1035.902 1057.623
## 93 1053.934 1043.703 1064.165 1038.287 1069.581
## 94 1060.947 1047.924 1073.969 1041.030 1080.863
## 95 1067.914 1052.301 1083.527 1044.036 1091.792
## 96 1074.867 1056.781 1092.953 1047.207 1102.527
## 97 1081.817 1061.325 1102.308 1050.478 1113.155
## 98 1088.765 1065.905 1111.625 1053.804 1123.726
## 99 1095.713 1070.501 1120.925 1057.154 1134.272
## 100 1102.661 1075.100 1130.222 1060.510 1144.812# Graficamos los pronósticos con los intervalos de predicción
autoplot(forecast_arima) +
labs(title = "Pronósticos ARIMA(1,2,2) para t=91 al t=100", x = "Tiempo", y = "yt") +
theme_minimal()