Taller 4 Topicos Cristian Suarez
2025-09-08
Ejercicio 1
# Cargar librerías
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries)
library(ggplot2)
library(readxl)
Mercado_laboral_y_población <- read_excel("Mercado laboral y población.xlsx")
head(Mercado_laboral_y_población)## # A tibble: 6 × 3
## Fecha `Tasa de desempleo - total nacional` Tasa de ocupación - total na…¹
## <chr> <chr> <chr>
## 1 yyyy/mm/dd % %
## 2 2001/01/31 16.62 57.58
## 3 2001/02/28 17.43 56.93
## 4 2001/03/31 15.81 57.57
## 5 2001/04/30 14.52 55.76
## 6 2001/05/31 14.04 56.23
## # ℹ abbreviated name: ¹`Tasa de ocupación - total nacional`datos <- Mercado_laboral_y_población
# 2. LIMPIAR Y CONVERTIR A NUMÉRICO
# Verificar el tipo de datos
cat("Tipo de datos de la columna:", class(datos$`Tasa de desempleo - total nacional`), "\n")## Tipo de datos de la columna: character# Convertir explícitamente a numérico y eliminar NA
valores_desempleo <- as.numeric(as.character(datos$`Tasa de desempleo - total nacional`))## Warning: NAs introduced by coercionvalores_desempleo <- valores_desempleo[!is.na(valores_desempleo)]
cat("Número de valores NA:", sum(is.na(valores_desempleo)), "\n")## Número de valores NA: 0cat("Longitud de la serie limpia:", length(valores_desempleo), "\n")## Longitud de la serie limpia: 295# 3. CREAR SERIE TEMPORAL NUMÉRICA
desempleo_ts <- ts(valores_desempleo,
start = c(2001, 1),
frequency = 12)
# a) GRÁFICO DE LA SERIE
plot(desempleo_ts, main = "Tasa de Desempleo Colombia (2001-2025)",
xlab = "Año", ylab = "Tasa de Desempleo (%)", col = "blue", lwd = 2)
grid()
# b) PRUEBA DE ESTACIONARIEDAD (ADF)
adf_test <- adf.test(desempleo_ts)
cat("\nPRUEBA ADF - Tasa de Desempleo:\n")##
## PRUEBA ADF - Tasa de Desempleo:cat("Dickey-Fuller =", round(adf_test$statistic, 4), "\n")## Dickey-Fuller = -3.1324cat("p-value =", adf_test$p.value, "\n")## p-value = 0.09959869if(adf_test$p.value > 0.05) {
cat("CONCLUSIÓN: La serie NO es estacionaria (p > 0.05)\n")
cat("Se requiere diferenciación para lograr estacionariedad\n")
} else {
cat("CONCLUSIÓN: La serie es estacionaria (p < 0.05)\n")
}## CONCLUSIÓN: La serie NO es estacionaria (p > 0.05)
## Se requiere diferenciación para lograr estacionariedad# c) DIFERENCIACIÓN - VERSIÓN SEGURA
desempleo_diff <- diff(desempleo_ts, differences = 1)
# Si aún hay error, usar esta alternativa:
# desempleo_diff <- na.omit(diff(desempleo_ts, differences = 1))
# d) IDENTIFICACIÓN DEL MODELO - FAC y FACP
par(mfrow = c(1, 2))
acf(desempleo_diff, na.action = na.pass, main = "FAC - Serie Diferenciada", lag.max = 24)
pacf(desempleo_diff, na.action = na.pass, main = "FACP - Serie Diferenciada", lag.max = 24)
par(mfrow = c(1, 1))
# e) MODELO SARIMA AUTOMÁTICO
cat("\nAJUSTANDO MODELO SARIMA...\n")##
## AJUSTANDO MODELO SARIMA...modelo <- auto.arima(desempleo_ts, seasonal = TRUE, stepwise = FALSE)
cat("\nMODELO SELECCIONADO:\n")##
## MODELO SELECCIONADO:print(modelo)## Series: desempleo_ts
## ARIMA(1,0,0)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 sma1 drift
## 0.9211 -0.8870 -0.0137
## s.e. 0.0237 0.0517 0.0086
##
## sigma^2 = 0.7017: log likelihood = -359.73
## AIC=727.46 AICc=727.6 BIC=742.04cat("\nRESUMEN DEL MODELO:\n")##
## RESUMEN DEL MODELO:summary(modelo)## Series: desempleo_ts
## ARIMA(1,0,0)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 sma1 drift
## 0.9211 -0.8870 -0.0137
## s.e. 0.0237 0.0517 0.0086
##
## sigma^2 = 0.7017: log likelihood = -359.73
## AIC=727.46 AICc=727.6 BIC=742.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01746578 0.8161182 0.5289384 -0.4196927 4.45886 0.420831
## ACF1
## Training set -0.10771# f) PRONÓSTICO (12 meses: Ago 2025 - Jul 2026)
pronostico <- forecast(modelo, h = 12)
# Gráfico del pronóstico
plot(pronostico, main = "Pronóstico Tasa de Desempleo (Ago 2025 - Jul 2026)",
xlab = "Año", ylab = "Tasa de Desempleo (%)")
lines(desempleo_ts, col = "black", lwd = 2)
grid()
# g) TABLA DE PRONÓSTICOS
cat("\nPRONÓSTICO TASA DE DESEMPLEO (%):\n")##
## PRONÓSTICO TASA DE DESEMPLEO (%):meses <- c("Ago-2025", "Sep-2025", "Oct-2025", "Nov-2025", "Dic-2025",
"Ene-2026", "Feb-2026", "Mar-2026", "Abr-2026", "May-2026",
"Jun-2026", "Jul-2026")
tabla_pronostico <- data.frame(
Mes = meses,
Pronostico = round(pronostico$mean, 2),
Limite_Inferior_80 = round(pronostico$lower[,1], 2),
Limite_Superior_80 = round(pronostico$upper[,1], 2),
Limite_Inferior_95 = round(pronostico$lower[,2], 2),
Limite_Superior_95 = round(pronostico$upper[,2], 2)
)
print(tabla_pronostico)## Mes Pronostico Limite_Inferior_80 Limite_Superior_80 Limite_Inferior_95
## 1 Ago-2025 8.34 7.27 9.42 6.70
## 2 Sep-2025 8.11 6.65 9.57 5.88
## 3 Oct-2025 7.73 6.00 9.45 5.09
## 4 Nov-2025 7.42 5.50 9.33 4.49
## 5 Dic-2025 8.04 5.97 10.10 4.88
## 6 Ene-2026 11.16 8.97 13.34 7.82
## 7 Feb-2026 9.99 7.71 12.27 6.51
## 8 Mar-2026 9.16 6.80 11.52 5.55
## 9 Abr-2026 9.50 7.08 11.92 5.79
## 10 May-2026 9.53 7.05 12.00 5.74
## 11 Jun-2026 9.14 6.62 11.66 5.28
## 12 Jul-2026 9.48 6.92 12.04 5.56
## Limite_Superior_95
## 1 9.99
## 2 10.34
## 3 10.36
## 4 10.34
## 5 11.19
## 6 14.50
## 7 13.48
## 8 12.77
## 9 13.21
## 10 13.32
## 11 12.99
## 12 13.39# h) VALIDACIÓN DEL MODELO
cat("\nVALIDACIÓN DEL MODELO:\n")##
## VALIDACIÓN DEL MODELO:cat("Prueba de Ljung-Box para residuales:\n")## Prueba de Ljung-Box para residuales:checkresiduals(modelo)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(0,1,1)[12] with drift
## Q* = 33.439, df = 22, p-value = 0.05595
##
## Model df: 2. Total lags used: 24# i) ECUACIÓN DEL MODELO
cat("\nECUACIÓN DEL MODELO SARIMA:\n")##
## ECUACIÓN DEL MODELO SARIMA:orden <- arimaorder(modelo)
cat("SARIMA(", orden[1], ",", orden[2], ",", orden[3], ")(",
orden[4], ",", orden[5], ",", orden[6], ")[12]\n", sep = "")## SARIMA(1,0,0)(0,1,1)[12]Ejercicio 3
# Cargar librerías
library(forecast)
library(tseries)
library(ggplot2)
library(readxl)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union# Cargar datos
X3_BEER <- read_excel("3. BEER.xlsx")
head(X3_BEER)## # A tibble: 6 × 2
## Quarter beer
## <chr> <dbl>
## 1 1956Q1 284.
## 2 1956Q2 213.
## 3 1956Q3 227.
## 4 1956Q4 308.
## 5 1957Q1 262
## 6 1957Q2 228.beer_data <- X3_BEER
# Convertir la columna Quarter a formato de fecha
# Extraer año y trimestre
beer_data$Year <- as.numeric(substr(beer_data$Quarter, 1, 4))
beer_data$Qtr <- as.numeric(substr(beer_data$Quarter, 6, 6))
# Crear serie de tiempo (frecuencia = 4 para datos trimestrales)
beer_ts <- ts(beer_data$beer, start = c(1956, 1), frequency = 4)
# =============================================================================
# a) Gráfico de la serie de tiempo y descripción
# =============================================================================
plot(beer_ts, main = "Producción de Cerveza en Australia (1956-1994)",
xlab = "Año", ylab = "Producción (megalitros)", col = "darkred", lwd = 2)
grid()
#
# -----------------------------------------------------------------------------
# La serie muestra una clara TENDENCIA CRECIENTE a lo largo del tiempo y un
# PATRÓN ESTACIONAL MARCADO con picos consistentes en el cuarto trimestre de
# cada año. La producción ha aumentado desde aproximadamente 200 megalitros
# en 1956 hasta alrededor de 600 megalitros en 1994, con fluctuaciones
# estacionales predecibles.
# =============================================================================
# b) Evaluación de estacionalidad (seasonal graph)
# =============================================================================
# Crear gráfico de cajas por trimestre
beer_data$Quarter_Factor <- factor(beer_data$Qtr, labels = c("Q1", "Q2", "Q3", "Q4"))
boxplot(beer ~ Quarter_Factor, data = beer_data,
main = "Producción por Trimestre",
xlab = "Trimestre", ylab = "Producción (megalitros)",
col = "lightblue")
# -----------------------------------------------------------------------------
# El gráfico de cajas por trimestre confirma un FUERTE PATRÓN ESTACIONAL.
# El cuarto trimestre (Q4) consistentemente presenta la producción más alta,
# seguido por el primer trimestre (Q1), mientras que el segundo y tercer
# trimestre (Q2 y Q3) muestran los niveles más bajos de producción.
# Esto refleja un comportamiento típico de consumo de cerveza con aumentos
# en periodos festivos.
# =============================================================================
# c) Prueba de estacionariedad (ADF)
# =============================================================================
adf_test <- adf.test(beer_ts)## Warning in adf.test(beer_ts): p-value greater than printed p-valueprint(adf_test)##
## Augmented Dickey-Fuller Test
##
## data: beer_ts
## Dickey-Fuller = 0.57302, Lag order = 5, p-value = 0.99
## alternative hypothesis: stationary# Mostrar resultados de la prueba
if(adf_test$p.value > 0.05) {
cat("La serie NO es estacionaria (p-value =", round(adf_test$p.value, 4), ")\n")
} else {
cat("La serie es estacionaria (p-value =", round(adf_test$p.value, 4), ")\n")
}## La serie NO es estacionaria (p-value = 0.99 )# -----------------------------------------------------------------------------
# Resultado: p-value = 0.99
# Conclusión: La serie NO ES ESTACIONARIA (p-value > 0.05).
# La alta probabilidad indica que no podemos rechazar la hipótesis nula de
# que existe una raíz unitaria, confirmando la presencia de tendencia y
# estacionalidad en la serie.
# =============================================================================
# d) Identificación y estimación del modelo SARIMA
# =============================================================================
# Diferenciación estacional
beer_diff <- diff(beer_ts, differences = 1)
beer_seasonal_diff <- diff(beer_ts, lag = 4, differences = 1)
# FAC y FACP de la serie diferenciada
par(mfrow = c(1, 2))
acf(beer_seasonal_diff, na.action = na.pass, main = "FAC - Serie diferenciada")
pacf(beer_seasonal_diff, na.action = na.pass, main = "FACP - Serie diferenciada")
par(mfrow = c(1, 1))
# Ajustar modelo SARIMA automáticamente
model <- auto.arima(beer_ts, seasonal = TRUE, stepwise = FALSE, approximation = FALSE)
summary(model)## Series: beer_ts
## ARIMA(0,1,2)(0,1,1)[4]
##
## Coefficients:
## ma1 ma2 sma1
## -0.9852 0.4156 -0.7668
## s.e. 0.0793 0.0814 0.0548
##
## sigma^2 = 256.8: log likelihood = -625.32
## AIC=1258.64 AICc=1258.92 BIC=1270.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4291945 15.60433 11.74947 -0.09639846 2.818065 0.7130395
## ACF1
## Training set 0.02234089# Mostrar ecuación del modelo
cat("\nEcuación del modelo SARIMA seleccionado:\n")##
## Ecuación del modelo SARIMA seleccionado:cat("SARIMA", arimaorder(model), "\n")## SARIMA 0 1 2 0 1 1 4# Modelo estimado: ARIMA(0,1,2)(0,1,1)[4]
# Ecuación del modelo:
# (1 - B)(1 - B⁴)Yₜ = (1 - 0.9852B + 0.4156B²)(1 - 0.7668B⁴)εₜ
# Justificación del modelo:
# - DIFERENCIACIÓN REGULAR (d=1): Necesaria para eliminar la tendencia
# - DIFERENCIACIÓN ESTACIONAL (D=1): Necesaria para eliminar la estacionalidad anual
# - TÉRMINOS MA(2): Los coeficientes MA son significativos (ma1 = -0.9852, ma2 = 0.4156)
# - TÉRMINO SMA(1): El coeficiente estacional es significativo (sma1 = -0.7668)
# - BONDAD DE AJUSTE: AIC = 1258.64, BIC = 1270.66 (valores relativamente bajos)
# - PRECISIÓN: MAPE = 2.82% (error porcentual medio bajo)
# =============================================================================
# e) Pronóstico hasta Q4 1996 (10 periodos hacia adelante)
# =============================================================================
forecast_beer <- forecast(model, h = 10)
# Gráfica del pronóstico
plot(forecast_beer, main = "Pronóstico de Producción de Cerveza hasta Q4 1996",
xlab = "Año", ylab = "Producción (megalitros)")
lines(beer_ts, col = "black", lwd = 2)
grid()
# Tabla de pronósticos
cat("\nPronóstico para 1994Q3 - 1996Q4:\n")##
## Pronóstico para 1994Q3 - 1996Q4:forecast_values <- data.frame(
Periodo = c("1994Q3", "1994Q4", "1995Q1", "1995Q2", "1995Q3", "1995Q4",
"1996Q1", "1996Q2", "1996Q3", "1996Q4"),
Pronostico = round(forecast_beer$mean, 2),
Lo_80 = round(forecast_beer$lower[,1], 2),
Hi_80 = round(forecast_beer$upper[,1], 2),
Lo_95 = round(forecast_beer$lower[,2], 2),
Hi_95 = round(forecast_beer$upper[,2], 2)
)
print(forecast_values)## Periodo Pronostico Lo_80 Hi_80 Lo_95 Hi_95
## 1 1994Q3 402.64 382.10 423.18 371.23 434.05
## 2 1994Q4 518.88 498.34 539.42 487.47 550.29
## 3 1995Q1 425.15 402.79 447.52 390.95 459.35
## 4 1995Q2 383.43 359.38 407.47 346.65 420.20
## 5 1995Q3 391.47 363.83 419.11 349.20 433.74
## 6 1995Q4 511.63 482.59 540.67 467.22 556.05
## 7 1996Q1 417.91 386.89 448.93 370.46 465.35
## 8 1996Q2 376.18 343.30 409.06 325.89 426.47
## 9 1996Q3 384.22 347.79 420.66 328.51 439.94
## 10 1996Q4 504.38 466.33 542.43 446.19 562.58# Pronósticos principales:
# - 1996Q4: 504.38 megalitros (rango: 446.19 - 562.58)
# - Patrón estacional mantenido: Los Q4 siguen siendo los trimestres con mayor producción
# Calidad del pronóstico: Los intervalos de confianza son relativamente estrechos
# (≈ ±40 megalitros para el 80% de confianza), indicando buena precisión predictiva.
# =============================================================================
# f) Prueba de residuales (Diagnóstico del modelo)
# =============================================================================
cat("\nPrueba de Ljung-Box para residuales:\n")##
## Prueba de Ljung-Box para residuales:checkresiduals(model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,1)[4]
## Q* = 3.2579, df = 5, p-value = 0.6603
##
## Model df: 3. Total lags used: 8# Prueba Ljung-Box: p-value = 0.6603
# Conclusión: Los residuales NO PRESENTAN AUTOCORRELACIÓN (p-value > 0.05),
# lo que confirma que el modelo captura adecuadamente la estructura de la serie
# y que los residuales son ruido blanco.
# El modelo SARIMA(0,1,2)(0,1,1)[4] es apropiado para la serie de producción
# de cerveza, capturando adecuadamente tanto la tendencia como la estacionalidad.
# Los pronósticos para 1995-1996 mantienen el patrón estacional histórico y
# presentan una precisión aceptable para la planificación de producción.