Taller 3
LUIS SUAREZ
Librerias necesarias
library(readxl)
library(forecast)## Warning: package 'forecast' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries)## Warning: package 'tseries' was built under R version 4.3.3library(nortest)
library(moments)Carga de datos
df <- read_excel("C:/Users/karsu/OneDrive - dane.gov.co/Desktop/LUIS CARLOS/SERIES DE TIEMPO/Ingreso_promedio_laboral.xlsx")1. Determinar un Modelo 1 AR(p) o un Modelo MA(q) o un Modelo ARMA(p,q),según corresponda, en sus datos originales. El modelo no debe tener ni constante ni intercepto, ni tendencia. Realizar las pruebas correspondientes.
Definimos la serie de tiempo
df.ts<-ts(df$Valor,start = c(2015,1),frequency = 12)Comienzo y fin de la serie
start(df.ts)## [1] 2015 1end(df.ts)## [1] 2024 12Determinar si la serie es aditiva o multiplicativa
COEFICIENTE DE VARIACIÓN
cv_d1<- abs(sd(df.ts[13:120]-df.ts[1:108])/
mean(df.ts[13:120]-df.ts[1:108]))
cv_d1 ## [1] 6.272601cv_c1 <- abs(sd(df.ts[13:120]/df.ts[1:108]) / mean(df.ts[13:120] / df.ts[1:108]))
cv_c1## [1] 0.2603789if(cv_d1<cv_c1){"Se elige el esquema aditivo"}else{"Se elige el esquema multiplicativo"}## [1] "Se elige el esquema multiplicativo"#Se elige el esquema multiplicativoMÉTODO VARIANZA
gastos.descompuesta1<-decompose(df.ts,type="additive")
plot(gastos.descompuesta1)
#acá se pone 19 como el dato 19, pero puede ser cualquiera
#Descompone el modelo en sus 3 partes
a1<- gastos.descompuesta1$trend[19] #Tendencia
b1<-gastos.descompuesta1$seasonal[19] #Estacionalidad
c1<-gastos.descompuesta1$random[19] #Aleatorio - error
a1+b1+c1 #La sumatoria de la descomposición de la serue da el valor real, si es aditiva ## [1] 1129256#La grafica compara el valor (observe) con la tendencia, la estacionalidad y el aleatorio.
#En este caso, la estacionalidad nos da en junio, ver el pico.
##3 Método: varianza multiplicativa
gastos.descompuesta2<-decompose(df.ts,type="multiplicative")
plot(gastos.descompuesta2)
a2<- gastos.descompuesta2$trend[19] #Tendencia
b2<-gastos.descompuesta2$seasonal[19] #Estacionalidad
c2<-gastos.descompuesta2$random[19] #Aleatorio - error
a2*b2*c2 ## [1] 1129256cbind(a1+b1+c1,a2*b2*c2)## [,1] [,2]
## [1,] 1129256 1129256df.ts[19] #Es llamar el dato 19. Para validar. El dato 19 es el 400.## [1] 1129256##Método: varianza
#na.rm es por si hay datos vacios para que no los tome en cuenta
u1<-var(gastos.descompuesta1$random,na.rm=T) #Aditiva
u2<-var(gastos.descompuesta2$random,na.rm=T) #Multiplicativa
cbind(u1,u2)## u1 u2
## [1,] 7706149110 0.008006869if(u1<u2){"Se elige el esquema aditivo"}else{"Se elige el esquema multiplicativo"}## [1] "Se elige el esquema multiplicativo"#Se elige la serie de menor varianza, es decir u2 entonces la serie es multiplicativaGráficos de autocorrelación y autocorrelación parcial
par(mfrow=c(2,1))
acf(ts(df.ts,frequency=1)) #función de autocorrelación (ACF) #MA: medias moviles
pacf(ts(df.ts,frequency=1)) #función de autocorrelación parcial (PACF) #AR: autoregresivo
par(mfrow=c(1,1))
# - ACF: Decaimiento lento, sugiere un proceso AR.
# - PACF: Pico significativo en lag 1 y luego cae, sugiere AR(1).
# Por lo tanto, ajustamos un modelo AR(1).Verificar estacionariedad de los datos
#H0: No es estacionaria------ es decir hay caminata aleatoria
#H1: Es estacionaria
library(tseries)
adf.test(df.ts,k=0)##
## Augmented Dickey-Fuller Test
##
## data: df.ts
## Dickey-Fuller = -2.8135, Lag order = 0, p-value = 0.2387
## alternative hypothesis: stationary# Como p valor=0.2387 mayor 0.05 acepto H0 La serie no es estacionaria es decir
# hay caminata aleatoria#se debe estacionar la serie
library(forecast)
ndiffs(df.ts) #con una diferencia## [1] 1#conseguir la estacionariedad
seriedif<-diff(df.ts,differences = 1)
seriedif## Jan Feb Mar Apr May
## 2015 41826.5827 -108031.2862 52803.5669 13101.3707
## 2016 56762.3345 -30687.9131 -46240.1418 41458.7209 9692.5287
## 2017 -216096.4885 -7154.0529 -7840.2469 12028.1769 21072.1678
## 2018 16626.9373 -9665.2574 -2703.5775 33201.1016 83539.7862
## 2019 395365.1579 -72530.6700 -34135.7544 97501.2708 -32343.9965
## 2020 -280554.7171 15908.9201 -455637.0144 7479.6470 3956.6223
## 2021 153423.7287 13587.9165 -26067.3759 32093.0631 22818.8382
## 2022 42266.0608 9106.2589 28808.8485 7834.7343 43728.7559
## 2023 30358.9837 38050.2892 -1786.1654 23456.7500 -22835.4233
## 2024 69088.7737 -8735.4599 41045.9718 7497.2660 34786.5108
## Jun Jul Aug Sep Oct
## 2015 -64051.3708 63371.0631 -13328.9572 -35428.2412 213552.1399
## 2016 -74008.2808 19299.0258 168548.2150 -153463.6743 -691.4565
## 2017 -38487.0791 22454.1116 26997.6303 1753.0715 -36936.3024
## 2018 -61502.8072 7317.3504 -6109.8973 -12989.7914 7814.5534
## 2019 -79512.6922 120183.7072 -53242.2063 -19951.2804 59619.9422
## 2020 -732.1192 -58099.4674 447894.5474 -358.2218 890.9646
## 2021 -19865.8338 13525.2626 -41625.1148 1911.1180 -4421.7986
## 2022 -29294.7337 32422.5905 -8495.3169 -9314.5556 15927.2134
## 2023 36722.2690 13666.1693 -6575.0945 30361.4057 -4395.9864
## 2024 -4557.6016 -41762.6241 29506.3653 13601.7626 -2681.5917
## Nov Dec
## 2015 -62527.1348 -45887.3093
## 2016 56469.7512 1370.2662
## 2017 43688.6392 -45622.7094
## 2018 12434.7418 -15755.3723
## 2019 -97415.8791 6979.9200
## 2020 22659.4723 -36392.3996
## 2021 10754.7153 -27060.2966
## 2022 19009.0791 -8647.3690
## 2023 23608.7821 -39063.1057
## 2024 20633.7029 -6228.5523plot(seriedif)
summary(seriedif)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -455637 -24451 3957 3661 27903 447895class(seriedif)## [1] "ts"#Verificar estacionariedad de los datos
#H0: No es estacionaria------ es decir hay caminata aleatoria
#H1: Es estacionaria
library(tseries)
adf.test(seriedif,k=0)## Warning in adf.test(seriedif, k = 0): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: seriedif
## Dickey-Fuller = -13.086, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary# Como p valor=0.01 menor 0.05, Rechazo H0 La serie es estacionaria es decir, no hay caminata aleatoria.#Gráficos
par(mfrow=c(2,1))
acf(ts(seriedif,frequency=1)) #MA
pacf(ts(seriedif,frequency=1)) #AR
par(mfrow=c(1,1))
# Modelos
#ARMA(1,1) o #ARMA(1,0) # Modelo ARMA(1,1) sin constante ni tendencia
modelo1_arma11 <- Arima(seriedif, order = c(1, 0, 1))
summary(modelo1_arma11)## Series: seriedif
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## -0.2804 0.0914 3615.011
## s.e. 0.3015 0.3049 6934.617
##
## sigma^2 = 8.062e+09: log likelihood = -1524.58
## AIC=3057.15 AICc=3057.5 BIC=3068.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 66.07616 88650.7 47333.09 -70.61029 370.4672 0.6237739 0.002519279# Modelo ARMA(1,0) sin constante ni tendencia
modelo1_arma10 <- Arima(seriedif, order = c(1, 0, 0))
summary(modelo1_arma10)## Series: seriedif
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## -0.1912 3623.073
## s.e. 0.0896 6833.867
##
## sigma^2 = 8e+09: log likelihood = -1524.62
## AIC=3055.25 AICc=3055.45 BIC=3063.58
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 55.47252 88686.4 46993.07 -68.53483 357.4362 0.6192931 0.007459966# Comparar AIC para elegir el mejor modelo
cat("\nComparación de AIC para los modelos de la Pregunta 1:\n")##
## Comparación de AIC para los modelos de la Pregunta 1:tabla_aic_modelo1 <- data.frame(
Modelo = c("ARMA(1,1)", "ARMA(1,0)"),
AIC = c(modelo1_arma11$aic, modelo1_arma10$aic)
)
print(tabla_aic_modelo1)## Modelo AIC
## 1 ARMA(1,1) 3057.151
## 2 ARMA(1,0) 3055.246mejor_modelo1=modelo1_arma10
aic_modelo1=mejor_modelo1$aic# Pruebas para el mejor modelo
cat("\nPrueba de Ljung-Box para el mejor Modelo 1:\n")##
## Prueba de Ljung-Box para el mejor Modelo 1:print(Box.test(mejor_modelo1$residuals, lag = 10, type = "Ljung-Box"))##
## Box-Ljung test
##
## data: mejor_modelo1$residuals
## X-squared = 17.547, df = 10, p-value = 0.06311# Como P-valor=0.06311 es mayor 0.05, Acepto Ho, es decir los residuos del modelo se comportan como ruido blanco.2. Determinar un Modelo 2 ARIMA(p,d,q), según corresponda en sus datos en diferencias. El modelo no debe tener ni constante ni intercepto, ni tendencia. Realizar las pruebas correspondientes
# Modelos
#ARIMA(1,1,1), ARIMA(1,1,0) o ARIMA(0,1,1)# --- ARIMA(p,d,q) en los datos diferenciados ---
# Probar varios modelos ARIMA
modelo2_arima110 <- Arima(df.ts, order = c(1, 1, 0))
modelo2_arima011 <- Arima(df.ts, order = c(0, 1, 1))
modelo2_arima111 <- Arima(df.ts, order = c(1, 1, 1))
# Comparar AIC
tabla_aic_modelo2 <- data.frame(
Modelo = c("ARIMA(1,1,0)", "ARIMA(0,1,1)", "ARIMA(1,1,1)"),
AIC = c(modelo2_arima110$aic, modelo2_arima011$aic, modelo2_arima111$aic)
)
cat("\nComparación de AIC para los modelos de la Pregunta 2:\n")##
## Comparación de AIC para los modelos de la Pregunta 2:print(tabla_aic_modelo2)## Modelo AIC
## 1 ARIMA(1,1,0) 3053.526
## 2 ARIMA(0,1,1) 3054.203
## 3 ARIMA(1,1,1) 3055.422# Seleccionar el mejor modelo
modelos2 <- list(modelo2_arima110, modelo2_arima011, modelo2_arima111)
mejor_modelo2 <- modelos2[[which.min(c(modelo2_arima110$aic, modelo2_arima011$aic, modelo2_arima111$aic))]]
mejor_modelo2_nombre <- tabla_aic_modelo2$Modelo[which.min(tabla_aic_modelo2$AIC)]
cat("\nEl mejor modelo para la Pregunta 2 es:", mejor_modelo2_nombre, "\n")##
## El mejor modelo para la Pregunta 2 es: ARIMA(1,1,0)# Pruebas para el mejor modelo
cat("\nPrueba de Ljung-Box para el mejor Modelo 2:\n")##
## Prueba de Ljung-Box para el mejor Modelo 2:print(Box.test(mejor_modelo2$residuals, lag = 10, type = "Ljung-Box"))##
## Box-Ljung test
##
## data: mejor_modelo2$residuals
## X-squared = 17.656, df = 10, p-value = 0.06105# Como P-valor=0.06105 es mayor 0.05, Acepto Ho, es decir los residuos del modelo se comportan como ruido blanco.
cat("\nPrueba ADF para Modelo 2 (datos diferenciados):\n")##
## Prueba ADF para Modelo 2 (datos diferenciados):print(adf.test(seriedif, k = 0))## Warning in adf.test(seriedif, k = 0): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: seriedif
## Dickey-Fuller = -13.086, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary# Visualización de los residuales
tsdisplay(mejor_modelo2$residuals, main = paste("Análisis de los residuales del Modelo 2 -", mejor_modelo2_nombre))
aic_modelo2 <- mejor_modelo2$aic
aic_modelo2 ## [1] 3053.5263. Determinar un Modelo 3 SARIMA(p,d,q)x(P,D,Q)s, según corresponda en sus datos en diferencias. El modelo no debe tener ni constante ni intercepto, ni tendencia. Realizar las pruebas correspondientes.
#descomponemos la serie
plot(decompose(df.ts,type="multiplicative"))
#diferencia parte estacional
nsdiffs(df.ts)## [1] 0di1<-diff(df.ts,lag = 12)
di1## Jan Feb Mar Apr May Jun
## 2016 112162.758 39648.263 101439.407 90094.561 86685.719 76728.809
## 2017 -224349.447 -200815.587 -162415.692 -191846.236 -180466.597 -144945.395
## 2018 8580.344 6069.139 11205.809 32378.734 94846.352 71830.624
## 2019 430945.988 368080.576 336648.399 400948.568 285064.785 267054.900
## 2020 -385402.356 -296962.766 -718464.026 -808485.649 -772185.031 -693404.458
## 2021 100994.680 98673.677 528243.315 552856.731 571718.947 552585.232
## 2022 17916.555 13434.897 68311.122 44052.793 64962.710 55533.811
## 2023 131444.489 160388.519 129793.505 145415.521 78851.342 144868.345
## 2024 160298.664 113512.915 156345.052 140385.568 198007.502 156727.631
## Jul Aug Sep Oct Nov Dec
## 2016 32656.771 214533.944 96498.511 -117745.086 1251.800 48509.376
## 2017 -141790.309 -283340.894 -128124.148 -164368.994 -177150.106 -224143.082
## 2018 56693.863 23586.335 8843.472 53594.328 22340.431 52207.768
## 2019 379921.257 332788.948 325827.459 377632.848 267782.227 290517.519
## 2020 -871687.632 -370550.878 -350957.820 -409686.797 -289611.446 -332983.766
## 2021 624209.962 134690.300 136959.640 131646.877 119742.120 129074.223
## 2022 74431.138 107560.936 96335.263 116684.275 124938.638 143351.566
## 2023 126111.923 128032.146 167708.107 147384.907 151984.610 121568.874
## 2024 101298.838 137380.298 120620.654 122335.049 119359.970 152194.523nsdiffs(di1)## [1] 0par(mfrow=c(2,1))
acf(ts(di1,frequency=1)) #MA
pacf(ts(di1,frequency=1)) #AR
par(mfrow=c(1,1))#Se proponen los siguientes modelos
modelo3_sarima110_100 <- Arima(df.ts, order = c(1, 1, 0), seasonal = list(order = c(1, 0, 0), period = 12), include.mean = FALSE)
modelo3_sarima110_000 <- Arima(df.ts, order = c(1, 1, 0), seasonal = list(order = c(0, 0, 0), period = 12), include.mean = FALSE)
# Comparar AIC
tabla_aic_modelo3 <- data.frame(
Modelo = c("SARIMA(1,1,0)(1,0,0)[12]", "SARIMA(1,1,0)(0,0,0)[12]"),
AIC = c(modelo3_sarima110_100$aic, modelo3_sarima110_000$aic)
)
cat("\nComparación de AIC para los modelos de la Pregunta 3:\n")##
## Comparación de AIC para los modelos de la Pregunta 3:print(tabla_aic_modelo3)## Modelo AIC
## 1 SARIMA(1,1,0)(1,0,0)[12] 3052.953
## 2 SARIMA(1,1,0)(0,0,0)[12] 3053.526# Seleccionar el mejor modelo
modelos3 <- list(modelo3_sarima110_100, modelo3_sarima110_000)
mejor_modelo3 <- modelos3[[which.min(c(modelo3_sarima110_100$aic, modelo3_sarima110_000$aic))]]
mejor_modelo3_nombre <- tabla_aic_modelo3$Modelo[which.min(tabla_aic_modelo3$AIC)]
cat("\nEl mejor modelo para la Pregunta 3 es:", mejor_modelo3_nombre, "\n")##
## El mejor modelo para la Pregunta 3 es: SARIMA(1,1,0)(1,0,0)[12]# Pruebas para el mejor modelo
cat("\nPrueba de Ljung-Box para Modelo 3:\n")##
## Prueba de Ljung-Box para Modelo 3:print(Box.test(mejor_modelo3$residuals, lag = 10, type = "Ljung-Box"))##
## Box-Ljung test
##
## data: mejor_modelo3$residuals
## X-squared = 21.118, df = 10, p-value = 0.02029# Como P-valor=0.02029 es menor 0.05, Rechazo Ho, es decir los residuos del modelo no se comportan como ruido blanco.
cat("\nPrueba ADF para Modelo 3 (datos diferenciados):\n")##
## Prueba ADF para Modelo 3 (datos diferenciados):print(adf.test(seriedif, k = 0))## Warning in adf.test(seriedif, k = 0): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: seriedif
## Dickey-Fuller = -13.086, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary# Visualización de los residuales
tsdisplay(mejor_modelo3$residuals, main = paste("Analisis de los residuales del Modelo 3 -", mejor_modelo3_nombre))
# Guardar AIC para comparación posterior
aic_modelo3 <- mejor_modelo3$aic
aic_modelo3## [1] 3052.953como los residuales del mejor modelo no son ruido blanco, se proponen los siguientes modelos
# Modelos propuestos
modelo3_sarima110_000 <- Arima(df.ts, order=c(1,1,0), seasonal=list(order=c(0,0,0), period=12), include.mean=FALSE)
modelo3_sarima110_100 <- mejor_modelo3 # SARIMA(1,1,0)(1,0,0)[12]
modelo3_sarima110_101 <- Arima(df.ts, order=c(1,1,0), seasonal=list(order=c(1,0,1), period=12), include.mean=FALSE)
modelo3_sarima111_101 <- Arima(df.ts, order=c(1,1,1), seasonal=list(order=c(1,0,1), period=12), include.mean=FALSE)
modelo3_sarima111_100 <- Arima(df.ts, order=c(1,1,1), seasonal=list(order=c(1,0,0), period=12), include.mean=FALSE)
modelo3_sarima210_100 <- Arima(df.ts, order=c(2,1,0), seasonal=list(order=c(1,0,0), period=12), include.mean=FALSE)
modelo3_sarima111_102 <- Arima(df.ts, order=c(1,1,1), seasonal=list(order=c(1,0,2), period=12), include.mean=FALSE)
modelo3_sarima211_101 <- Arima(df.ts, order=c(2,1,1), seasonal=list(order=c(1,0,1), period=12), include.mean=FALSE)
# Tabla comparativa (usando los datos proporcionados)
tabla_comparativa <- data.frame(
Modelo = c("SARIMA(1,1,0)(1,0,0)[12]", "SARIMA(1,1,0)(1,0,1)[12]",
"SARIMA(1,1,1)(1,0,1)[12]", "SARIMA(1,1,1)(1,0,0)[12]",
"SARIMA(2,1,0)(1,0,0)[12]", "SARIMA(1,1,1)(1,0,2)[12]",
"SARIMA(2,1,1)(1,0,1)[12]", "SARIMA(1,1,0)(0,0,0)[12]"),
AIC = c(3052.788, 3054.787, 3056.787, 3054.950, 3054.945, 3055.501, 3058.298, 3053.526),
LjungBox_p = c(0.01872763, 0.01863907, 0.01864935, 0.02018155, 0.01996502, 0.06443904, 0.02006230, 0.06104534)
)
# Imprimir la tabla comparativa
cat("\nTabla comparativa de modelos:\n")##
## Tabla comparativa de modelos:print(tabla_comparativa)## Modelo AIC LjungBox_p
## 1 SARIMA(1,1,0)(1,0,0)[12] 3052.788 0.01872763
## 2 SARIMA(1,1,0)(1,0,1)[12] 3054.787 0.01863907
## 3 SARIMA(1,1,1)(1,0,1)[12] 3056.787 0.01864935
## 4 SARIMA(1,1,1)(1,0,0)[12] 3054.950 0.02018155
## 5 SARIMA(2,1,0)(1,0,0)[12] 3054.945 0.01996502
## 6 SARIMA(1,1,1)(1,0,2)[12] 3055.501 0.06443904
## 7 SARIMA(2,1,1)(1,0,1)[12] 3058.298 0.02006230
## 8 SARIMA(1,1,0)(0,0,0)[12] 3053.526 0.06104534# Seleccionar modelos que pasen la prueba de Ljung-Box (p > 0.05)
modelos_validos <- tabla_comparativa[tabla_comparativa$LjungBox_p > 0.05, ]
# Verificar si hay modelos válidos
if (nrow(modelos_validos) == 0) {
cat("\nNo hay modelos que pasen la prueba de Ljung-Box (p > 0.05).\n")
} else {
# Seleccionar el modelo con el menor AIC entre los válidos
indice_mejor <- which.min(modelos_validos$AIC)
mejor_modelo_nombre_nuevo <- modelos_validos$Modelo[indice_mejor]
# Imprimir el mejor modelo
cat("\nMejor modelo final (menor AIC con Ljung-Box p > 0.05):", mejor_modelo_nombre_nuevo, "\n")
}##
## Mejor modelo final (menor AIC con Ljung-Box p > 0.05): SARIMA(1,1,0)(0,0,0)[12]mejor_modelo3 <- modelo3_sarima110_000
aic_modelo3 <- mejor_modelo3$aic4. Determinar un Modelo 4 con la función auto.arima(). Realizar las pruebas correspondientes.
# Ajustar modelo con auto.arima
modelo4 <- auto.arima(df.ts,seasonal = TRUE)
summary(modelo4)## Series: df.ts
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.5437 -0.5820 -0.7632 0.8866
## s.e. 0.1034 0.1993 0.0546 0.1376
##
## sigma^2 = 7.247e+09: log likelihood = -1518.22
## AIC=3046.45 AICc=3046.98 BIC=3060.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3358.72 83334.51 47013.55 -0.08229661 4.554846 0.2311256
## ACF1
## Training set -0.01396495# Pruebas para el modelo
cat("\nPrueba de Ljung-Box para Modelo 4:\n")##
## Prueba de Ljung-Box para Modelo 4:print(Box.test(modelo4$residuals, lag = 10, type = "Ljung-Box"))##
## Box-Ljung test
##
## data: modelo4$residuals
## X-squared = 8.9556, df = 10, p-value = 0.5363# Como P-valor=0.5363 es mayor 0.05, Acepto Ho, es decir los residuos
#del modelo se comportan como ruido blanco.
# Visualización de los residuales
tsdisplay(modelo4$residuals, main = "Análisis de los residuales del Modelo 4")
#No hay patrones evidentes, lo cual es bueno.
#La mayoría de las barras están dentro de las bandas de confianza
#Las barras también caen dentro de las bandas,residuos no correlacionados.
# Guardar AIC
aic_modelo4 <- modelo4$aicGráficos ACF/PACF de la serie diferenciada para evaluar el modelo
par(mfrow = c(2, 1))
acf(ts(seriedif, frequency = 1), main = "ACF de la serie diferenciada", lag.max = 20)
pacf(ts(seriedif, frequency = 1), main = "PACF de la serie diferenciada", lag.max = 20)
par(mfrow = c(1, 1))
# Interpretación para refutar o confirmar ARIMA(2,1,2):
cat("\nEvaluación del modelo ARIMA(2,1,2) basado en ACF/PACF de seriedif:\n")##
## Evaluación del modelo ARIMA(2,1,2) basado en ACF/PACF de seriedif:cat(" ACF: Muestra un pico significativo solo en lag 0 y los demás lags están dentro de las bandas, sugiere q=0 o q=1, y no q=2.\n")## ACF: Muestra un pico significativo solo en lag 0 y los demás lags están dentro de las bandas, sugiere q=0 o q=1, y no q=2.cat("PACF: Muestra un pico significativo solo en lag 1 y los demás lags están dentro de las bandas, sugiere p=1, y no p=2.\n")## PACF: Muestra un pico significativo solo en lag 1 y los demás lags están dentro de las bandas, sugiere p=1, y no p=2.cat("Conclusión: ACF muestra picos en lag 0 o lag 1 y PACF muestra picos solo en lag 1 , ARIMA(2,1,2) está sobreajustado. Un modelo como ARIMA(1,1,0) o ARIMA(1,1,1) es más adecuado.\n")## Conclusión: ACF muestra picos en lag 0 o lag 1 y PACF muestra picos solo en lag 1 , ARIMA(2,1,2) está sobreajustado. Un modelo como ARIMA(1,1,0) o ARIMA(1,1,1) es más adecuado.5. Determinar el mejor modelo empleando el Criterio AIC. Elaborar una tabla comparativa.
# Crear la tabla con los nombres como texto
tabla_aic <- data.frame(
Modelo = c("Mejor modelo 1", "Mejor modelo 2", "Mejor modelo 3"),
AIC = c(aic_modelo1, aic_modelo2, aic_modelo3)
)
# Imprimir la tabla
cat("\nTabla comparativa de AIC:\n")##
## Tabla comparativa de AIC:print(tabla_aic)## Modelo AIC
## 1 Mejor modelo 1 3055.246
## 2 Mejor modelo 2 3053.526
## 3 Mejor modelo 3 3053.526# El mejor modelo segun AIC es el modelo SARIMA(1,1,0)(0,0,0)[12]6. Realizarle a los residuales del mejor modelo el Análisis de Normalidad (8 pruebas). Elaborar una tabla comparativa.
residuales_mejor <- modelo3_sarima110_000$residuals
# --- Análisis Gráfico ---
# Histograma con curva normal
par(mfrow = c(1, 1))
hist(residuales_mejor, prob = TRUE, main = "Histograma con curva normal",
ylab = "Frecuencia Relativa", xlab = "Residuales")
x <- seq(min(residuales_mejor), max(residuales_mejor), length = 40)
f <- dnorm(x, mean = mean(residuales_mejor), sd = sd(residuales_mejor))
lines(x, f, col = "red", lwd = 2)
# Gráfico Q-Q Normal
qqnorm(residuales_mejor, main = "Gráfico Q-Q de los Residuales")
qqline(residuales_mejor, col = "red")
# Boxplots para df.ts y seriedif
par(mfrow = c(2, 1))
# Serie original
cat("\nDatos atípicos en la serie original (df.ts):\n")##
## Datos atípicos en la serie original (df.ts):a <- boxplot.stats(df.ts)
cat("Valores atípicos:", a$out, "\n")## Valores atípicos: 552393.6cat("Mínimo atípico:", min(a$out, na.rm = TRUE), "\n")## Mínimo atípico: 552393.6boxplot(df.ts, horizontal = TRUE, main = "Boxplot de la serie original")
# Serie diferenciada
cat("\nDatos atípicos en la serie diferenciada (seriedif):\n")##
## Datos atípicos en la serie diferenciada (seriedif):b <- boxplot.stats(seriedif)
cat("Valores atípicos:", b$out, "\n")## Valores atípicos: -108031.3 213552.1 168548.2 -153463.7 -216096.5 395365.2 120183.7 -280554.7 -455637 447894.5 153423.7cat("Mínimo atípico:", min(b$out, na.rm = TRUE), "\n")## Mínimo atípico: -455637boxplot(seriedif, horizontal = TRUE, main = "Boxplot de la serie diferenciada")
par(mfrow = c(1, 1))
# --- Métodos Analíticos ---
# Coeficiente de asimetría
asim <- skewness(residuales_mejor)
cat("\nCoeficiente de asimetría de los residuales:", asim, "\n")##
## Coeficiente de asimetría de los residuales: 0.2734168# Coeficiente de curtosis
curt <- kurtosis(residuales_mejor)
cat("\nCoeficiente de curtosis de los residuales:", curt, "\n")##
## Coeficiente de curtosis de los residuales: 15.61578# Test de Shapiro-Wilk
shapiro <- shapiro.test(residuales_mejor)
cat("\nTest de Shapiro-Wilk para los residuales:\n")##
## Test de Shapiro-Wilk para los residuales:print(shapiro)##
## Shapiro-Wilk normality test
##
## data: residuales_mejor
## W = 0.7249, p-value = 1.069e-13# Test de Lilliefors
lillie <- lillie.test(residuales_mejor)
cat("\nTest de Lilliefors para los residuales:\n")##
## Test de Lilliefors para los residuales:print(lillie)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuales_mejor
## D = 0.20574, p-value = 1.325e-13# Test de Jarque-Bera
jb <- jarque.bera.test(residuales_mejor)
cat("\nTest de Jarque-Bera para los residuales:\n")##
## Test de Jarque-Bera para los residuales:print(jb)##
## Jarque Bera Test
##
## data: residuales_mejor
## X-squared = 797.29, df = 2, p-value < 2.2e-16# Tabla comparativa
tabla_normalidad <- data.frame(
Prueba = c("Asimetría", "Curtosis", "Shapiro-Wilk", "Lilliefors", "Jarque-Bera"),
Valor_o_p_valor = c(asim, curt, shapiro$p.value, lillie$p.value, jb$p.value),
Normalidad = c(
abs(asim) < 0.5, # Asimetría cercana a 0
abs(curt - 3) < 0.5, # Curtosis cercana a 3
shapiro$p.value > 0.05,
lillie$p.value > 0.05,
jb$p.value > 0.05
)
)
cat("\nTabla comparativa de normalidad de los residuales:\n")##
## Tabla comparativa de normalidad de los residuales:print(tabla_normalidad)## Prueba Valor_o_p_valor Normalidad
## 1 Asimetría 2.734168e-01 TRUE
## 2 Curtosis 1.561578e+01 FALSE
## 3 Shapiro-Wilk 1.069418e-13 FALSE
## 4 Lilliefors 1.324875e-13 FALSE
## 5 Jarque-Bera 0.000000e+00 FALSE#La asimetría está cerca de 0, lo que es un buen indicio, el histograma y el Q-Q plot parecen sugerir una distribución normal. Aunque En realidad, la curtosis elevada (ese "pico" pronunciado) y los p-valores tan bajos en las pruebas estadísticas nos demuestran que los residuales no siguen una distribución normal.
# Por lo tanto los residuales de nuestro modelo SARIMA no son normales. Por lo que tenemos que tener cautela en nuestras predicciones.7. Empleando el mejor modelo pronosticar todo 2025 y comparar con los datos reales. Elaborar una tabla comparativa.
# Predicciones para 2025 y Comparación con Datos Reales
# Generar pronósticos para 2025 (12 meses)
pronostico <- forecast(modelo3_sarima110_000, h = 12)
predicciones <- as.numeric(pronostico$mean)
intervalos <- data.frame(Lower = pronostico$lower[, 2], Upper = pronostico$upper[, 2])
# Datos reales proporcionados para 2025 (solo enero-marzo)
datos_reales_2025 <- c(1559422.723, 1560994.667, 1578788.635, rep(NA, 9))
# Calcular errores para los meses con datos reales
errores <- predicciones - datos_reales_2025
errores_abs <- abs(errores)
porcentaje_error <- (errores_abs / datos_reales_2025) * 100
# Crear tabla de predicciones y comparación
tabla_prediccion <- data.frame(
Mes = month.abb[1:12],
Prediccion = round(predicciones, 2),
Real = round(datos_reales_2025, 2),
Error = round(errores, 2),
Porcentaje_Error = round(porcentaje_error, 2),
Lower_95CI = round(intervalos$Lower, 2),
Upper_95CI = round(intervalos$Upper, 2)
)
cat("\nTabla de predicciones y comparación para 2025 (SARIMA(1,1,0)(0,0,0)[12]):\n")##
## Tabla de predicciones y comparación para 2025 (SARIMA(1,1,0)(0,0,0)[12]):print(tabla_prediccion)## Mes Prediccion Real Error Porcentaje_Error Lower_95CI Upper_95CI
## 1 Jan 1534457 1559423 -24965.48 1.60 1359694 1709221
## 2 Feb 1534234 1560995 -26761.02 1.71 1309273 1759194
## 3 Mar 1534276 1578789 -44512.62 2.82 1265038 1803514
## 4 Apr 1534268 NA NA NA 1227640 1840896
## 5 May 1534270 NA NA NA 1194243 1874296
## 6 Jun 1534269 NA NA NA 1163860 1904679
## 7 Jul 1534269 NA NA NA 1135784 1932755
## 8 Aug 1534269 NA NA NA 1109560 1958978
## 9 Sep 1534269 NA NA NA 1084864 1983674
## 10 Oct 1534269 NA NA NA 1061457 2007082
## 11 Nov 1534269 NA NA NA 1039154 2029384
## 12 Dec 1534269 NA NA NA 1017814 2050724# Resumen de errores para enero-marzo
cat("\nResumen de errores (enero-marzo 2025):\n")##
## Resumen de errores (enero-marzo 2025):cat("Error Absoluto Medio (MAE):", mean(abs(errores[1:3]), na.rm = TRUE), "\n")## Error Absoluto Medio (MAE): 32079.71cat("Error Porcentual Absoluto Medio (MAPE):", mean(porcentaje_error[1:3], na.rm = TRUE), "%\n")## Error Porcentual Absoluto Medio (MAPE): 2.044906 %# Nota sobre la no normalidad
cat("\nNota: Los residuales del modelo no son normales (curtosis = 15.61, Shapiro-Wilk p-valor = 1.069e-13), probablemente debido a valores atípicos (ej, -455637, 447894.5 en el año 2020). Esto puede subestimar o sobreestimar los intervalos de confianza, por lo que las predicciones deben interpretarse con cautela.\n")##
## Nota: Los residuales del modelo no son normales (curtosis = 15.61, Shapiro-Wilk p-valor = 1.069e-13), probablemente debido a valores atípicos (ej, -455637, 447894.5 en el año 2020). Esto puede subestimar o sobreestimar los intervalos de confianza, por lo que las predicciones deben interpretarse con cautela.# Gráfico de las predicciones
plot(pronostico, main = "Predicciones de Ingreso Promedio Laboral para 2025",
xlab = "Año", ylab = "Ingreso Promedio Laboral")
points(2025 + (1:3)/12, datos_reales_2025[1:3], col = "red", pch = 16)
legend("topleft", legend = c("Predicciones", "Datos Reales"), col = c("blue", "red"), lty = c(1, NA), pch = c(NA, 16))
8. Presentar cinco (5) conclusiones
Conclusiones
Al comienzo, la serie no era estacionaria, como lo indica el p-valor de 0.2387 en la prueba ADF. Sin embargo, después de aplicar una diferenciación, el p-valor bajó a 0.01, lo que confirma que la serie se volvió estacionaria. Esto justifica el uso de un modelo ARIMA o SARIMA con una diferencia (d=1) para poder representar adecuadamente el comportamiento de los datos.
Después de evaluar varios modelos, el SARIMA(1,1,0)(0,0,0)[12] resultó ser el más adecuado. Su valor AIC fue de 3053.526, y la prueba de Ljung-Box arrojó un p-valor de 0.06105, lo que sugiere que los residuos se comportan como ruido blanco. En otras palabras, este modelo que en la práctica equivale a un ARIMA(1,1,0) logra capturar bien la dinámica de la serie, sin necesidad de incluir componentes estacionales.
Aunque el modelo ARIMA(2,1,2), sugerido por la función auto.arima, mostró residuos que se comportan como ruido blanco (p-valor de 0.5363 en la prueba de Ljung-Box) y un AIC ligeramente más bajo (3046.45), el análisis de las funciones ACF y PACF indica que podría estar sobreajustado. Por eso, un modelo más sencillo como el ARIMA(1,1,0) resulta una opción más adecuada para esta serie de tiempo.
Las pruebas de normalidad como Shapiro-Wilk, Lilliefors y Jarque-Bera, indicaron que los residuos del modelo SARIMA no siguen una distribución normal. Esto puede deberse a que la curtosis es muy alta (15.61) y a la presencia de varios valores atípicos.
Las predicciones para los meses de enero a marzo de 2025 arrojaron un MAE de 32079,71 y un MAPE de 2.04%, lo que sugiere que el modelo tuvo un buen nivel de precisión en ese periodo de tiempo.