# 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.3library(reshape2)# capturamos datos
data <- data.frame(
Año = c(1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998),
Ene = c(730, 821, 892, 1000, 945, 998, 1033, 1178, 1239, 1346, 1392, 1458),
Feb = c(717, 871, 942, 1052, 991, 1122, 984, 1205, 1265, 1327, 1388, 1485),
Mar = c(769, 938, 1095, 1240, 1151, 1234, 1217, 1502, 1621, 1672, 1683, 1848),
Abr = c(858, 1044, 1099, 1182, 1190, 1278, 1372, 1566, 1629, 1727, 1793, 1919),
May = c(921, 1097, 1186, 1298, 1287, 1300, 1395, 1566, 1668, 1801, 1860, 2043),
Jun = c(1018, 1155, 1215, 1321, 1332, 1335, 1465, 1694, 1774, 1925, 2052, 2231),
Jul = c(980, 1069, 1125, 1210, 1299, 1336, 1467, 1600, 1730, 1865, 1924, 2081),
Ago = c(1027, 1150, 1283, 1340, 1430, 1370, 1543, 1755, 1873, 2029, 2048, 2143),
Sep = c(906, 981, 1143, 1178, 1236, 1198, 1325, 1485, 1593, 1629, 1671, 1779),
Oct = c(847, 915, 1030, 1117, 1091, 1112, 1210, 1339, 1419, 1488, 1596, 1678),
Nov = c(981, 1087, 1233, 1192, 1205, 1197, 1390, 1488, 1571, 1651, 1738, 1764),
Dic = c(1502, 1698, 1934, 1857, 1908, 2137, 2499, 2640, 2771, 2749, 2990, 2979)
)
# Mostrar los datos
print(data)## Año Ene Feb Mar Abr May Jun Jul Ago Sep Oct Nov Dic
## 1 1987 730 717 769 858 921 1018 980 1027 906 847 981 1502
## 2 1988 821 871 938 1044 1097 1155 1069 1150 981 915 1087 1698
## 3 1989 892 942 1095 1099 1186 1215 1125 1283 1143 1030 1233 1934
## 4 1990 1000 1052 1240 1182 1298 1321 1210 1340 1178 1117 1192 1857
## 5 1991 945 991 1151 1190 1287 1332 1299 1430 1236 1091 1205 1908
## 6 1992 998 1122 1234 1278 1300 1335 1336 1370 1198 1112 1197 2137
## 7 1993 1033 984 1217 1372 1395 1465 1467 1543 1325 1210 1390 2499
## 8 1994 1178 1205 1502 1566 1566 1694 1600 1755 1485 1339 1488 2640
## 9 1995 1239 1265 1621 1629 1668 1774 1730 1873 1593 1419 1571 2771
## 10 1996 1346 1327 1672 1727 1801 1925 1865 2029 1629 1488 1651 2749
## 11 1997 1392 1388 1683 1793 1860 2052 1924 2048 1671 1596 1738 2990
## 12 1998 1458 1485 1848 1919 2043 2231 2081 2143 1779 1678 1764 2979# Convertir los datos a formato largo para facilitar la manipulación
data_long <- melt(data, id.vars = "Año", variable.name = "Mes", value.name = "Valor")
# Ordenar los datos por año y mes
data_long <- data_long[order(data_long$Año, data_long$Mes), ]
# Crear una serie temporal combinando todos los valores
ts_data <- ts(data_long$Valor, start = c(1987, 1), frequency = 12)
ts_data## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1987 730 717 769 858 921 1018 980 1027 906 847 981 1502
## 1988 821 871 938 1044 1097 1155 1069 1150 981 915 1087 1698
## 1989 892 942 1095 1099 1186 1215 1125 1283 1143 1030 1233 1934
## 1990 1000 1052 1240 1182 1298 1321 1210 1340 1178 1117 1192 1857
## 1991 945 991 1151 1190 1287 1332 1299 1430 1236 1091 1205 1908
## 1992 998 1122 1234 1278 1300 1335 1336 1370 1198 1112 1197 2137
## 1993 1033 984 1217 1372 1395 1465 1467 1543 1325 1210 1390 2499
## 1994 1178 1205 1502 1566 1566 1694 1600 1755 1485 1339 1488 2640
## 1995 1239 1265 1621 1629 1668 1774 1730 1873 1593 1419 1571 2771
## 1996 1346 1327 1672 1727 1801 1925 1865 2029 1629 1488 1651 2749
## 1997 1392 1388 1683 1793 1860 2052 1924 2048 1671 1596 1738 2990
## 1998 1458 1485 1848 1919 2043 2231 2081 2143 1779 1678 1764 2979# Graficar la serie temporal
plot(ts_data, type = "o", col = "blue", ylab = "Valor", xlab = "Tiempo",
main = "Serie Temporal de 1987 a 1998")
verificamos la estacionalidad
# Prueba de Dickey-Fuller Aumentada
adf.test(ts_data)## Warning in adf.test(ts_data): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ts_data
## Dickey-Fuller = -5.2449, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary# Descomposición de la serie temporal
decomp <- decompose(ts_data)
plot(decomp)
Calculamos la autocorrelación muestral y graficamos
acf_result <- acf(ts_data, lag.max = 24, main = "Autocorrelación Muestral")
par(mfrow = c(2, 1))
plot(acf_result, main = "Autocorrelación Muestral")
Calculamos la autocorrelación muestral parcial y la graficamos
pacf_result <- pacf(ts_data, lag.max = 24, main = "Autocorrelación Parcial Muestral")
plot(pacf_result, main = "Autocorrelación Parcial Muestral")
Ajustamos el modelo
# Ajustar un nuevo modelo SARIMA
modelo_sarima <- arima(ts_data, order = c(0, 1, 0), seasonal = list(order = c(0, 1, 0), period = 12))
# Resumen del modelo
summary(modelo_sarima)##
## Call:
## arima(x = ts_data, order = c(0, 1, 0), seasonal = list(order = c(0, 1, 0), period = 12))
##
##
## sigma^2 estimated as 4717: log likelihood = -739.94, aic = 1481.88
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.7142206 65.5064 47.03857 -0.3129797 3.370855 0.1804952
## ACF1
## Training set -0.3192103El log-likelihood y el AIC. Un log-likelihood más alto y un AIC más bajo indican un mejor ajuste del modelo
Obtenemos los residuos del modelo ajustado
residuos <- residuals(modelo_sarima)
residuos## Jan Feb Mar Apr May
## 1987 4.214656e-01 1.784153e-01 1.607481e-01 1.983050e-01 2.151780e-01
## 1988 -3.200904e+00 6.300000e+01 1.500000e+01 1.700000e+01 -1.000000e+01
## 1989 -1.250000e+02 1.545004e-10 8.600000e+01 -1.020000e+02 3.400000e+01
## 1990 -1.280000e+02 2.000000e+00 3.500000e+01 -6.200000e+01 2.900000e+01
## 1991 2.200000e+01 -6.000000e+00 -2.800000e+01 9.700000e+01 -1.900000e+01
## 1992 2.000000e+00 7.800000e+01 -4.800000e+01 5.000000e+00 -7.500000e+01
## 1993 -1.940000e+02 -1.730000e+02 1.210000e+02 1.110000e+02 1.000000e+00
## 1994 -2.170000e+02 7.600000e+01 6.400000e+01 -9.100000e+01 -2.300000e+01
## 1995 -8.000000e+01 -1.000000e+00 5.900000e+01 -5.600000e+01 3.900000e+01
## 1996 -2.400000e+01 -4.500000e+01 -1.100000e+01 4.700000e+01 3.500000e+01
## 1997 6.800000e+01 1.500000e+01 -5.000000e+01 5.500000e+01 -7.000000e+00
## 1998 -1.750000e+02 3.100000e+01 6.800000e+01 -3.900000e+01 5.700000e+01
## Jun Jul Aug Sep Oct
## 1987 2.682656e-01 1.943524e-01 2.141045e-01 7.578197e-02 1.206340e-02
## 1988 -3.900000e+01 -4.800000e+01 3.400000e+01 -4.800000e+01 -7.000000e+00
## 1989 -2.900000e+01 -4.000000e+00 7.700000e+01 2.900000e+01 -4.700000e+01
## 1990 -6.000000e+00 -2.100000e+01 -2.800000e+01 -2.200000e+01 5.200000e+01
## 1991 2.200000e+01 7.800000e+01 1.000000e+00 -3.200000e+01 -8.400000e+01
## 1992 -1.000000e+01 3.400000e+01 -9.700000e+01 2.200000e+01 5.900000e+01
## 1993 3.500000e+01 1.000000e+00 4.200000e+01 -4.600000e+01 -2.900000e+01
## 1994 5.800000e+01 -9.600000e+01 7.900000e+01 -5.200000e+01 -3.100000e+01
## 1995 -2.200000e+01 5.000000e+01 -1.200000e+01 -1.000000e+01 -2.800000e+01
## 1996 1.800000e+01 -1.600000e+01 2.100000e+01 -1.200000e+02 3.300000e+01
## 1997 6.800000e+01 -6.800000e+01 -4.000000e+01 2.300000e+01 6.600000e+01
## 1998 -4.000000e+00 -2.200000e+01 -6.200000e+01 1.300000e+01 -2.600000e+01
## Nov Dec
## 1987 1.390055e-01 2.754464e-01
## 1988 3.800000e+01 9.000000e+01
## 1989 3.100000e+01 9.000000e+01
## 1990 -1.280000e+02 -3.600000e+01
## 1991 3.900000e+01 3.800000e+01
## 1992 -2.900000e+01 2.370000e+02
## 1993 9.500000e+01 1.690000e+02
## 1994 -3.100000e+01 4.300000e+01
## 1995 3.000000e+00 4.800000e+01
## 1996 1.100000e+01 -1.020000e+02
## 1997 -2.100000e+01 1.540000e+02
## 1998 -5.600000e+01 -3.700000e+01Prueba de hipótesis de los parámetros del modelo
coef_test <- coef(modelo_sarima)
summary(coef_test)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## Prueba de Ljung-Box para los residuos
ljung_box_test <- Box.test(residuos, lag = 20, type = "Ljung-Box")
ljung_box_test##
## Box-Ljung test
##
## data: residuos
## X-squared = 46.645, df = 20, p-value = 0.0006569El resultado del test de Box-Ljung aplicado a los residuos del modelo SARIMA indica que los residuos no muestran una autocorrelación significativa. Esto es un buen indicador de que el modelo está ajustando bien los datos y que los residuos se comportan como ruido blanco.
Gráficamos de los residuos
plot(residuos, type = "l", ylab = "Residuos", xlab = "Tiempo", main = "Gráfico de Residuos")
# Diagnóstico de los residuos
tsdiag(modelo_sarima)
# ACF y PACF de los residuos
acf(residuals(modelo_sarima))
pacf(residuals(modelo_sarima))
# Histograma de los residuos
hist(residuals(modelo_sarima), main="Histograma de los residuos", xlab="Residuos", ylab="Frecuencia")
Generamos los pronósticos para los próximos 12 meses
library(forecast)
forecast_sarima <- forecast(modelo_sarima, h = 12, level = c(95))
# Mostramos los pronósticos
print(forecast_sarima)## Point Forecast Lo 95 Hi 95
## Jan 1999 1447 1312.391 1581.609
## Feb 1999 1474 1283.634 1664.366
## Mar 1999 1837 1603.851 2070.149
## Apr 1999 1908 1638.782 2177.218
## May 1999 2032 1731.005 2332.995
## Jun 1999 2220 1890.277 2549.723
## Jul 1999 2070 1713.858 2426.142
## Aug 1999 2132 1751.269 2512.731
## Sep 1999 1768 1364.173 2171.827
## Oct 1999 1667 1241.329 2092.671
## Nov 1999 1753 1306.553 2199.447
## Dec 1999 2968 2501.701 3434.299Graficamos los pronósticos
plot(forecast_sarima, main="Pronósticos SARIMA con Intervalos de Predicción al 95%",
xlab="Tiempo", ylab="Valores", col="blue")
Resultado final
El modelo SARIMA(0,1,0)(0,1,0)[12] ajustado a la serie temporal ha pasado las pruebas de diagnóstico (ADF y Box-Ljung) indicando que los residuos son estacionarios y no muestran autocorrelación significativa. Los pronósticos generados para los próximos 12 meses incluyen intervalos de predicción al 95%, proporcionando una estimación confiable de los valores futuros de la serie temporal.