Módulo 7·Ingeniería, economía y ciencias ambientales

Ejemplos Completos de Aplicación

Tres casos de uso detallados con código R completo: demanda eléctrica horaria, PIB trimestral agrícola de Colombia y temperatura mensual de Nottingham Castle.

Demanda eléctrica horariaPIB trimestral (agricultura)Temperatura Nottingham CastleCódigo R reproducibleInterpretación de resultados

Anterior: Síntesis Metodológica Siguiente: Pruebas de Hipótesis y Gráficos Diagnósticos

¿Por qué ver ejemplos completos?

La teoría de los módulos 1–6 dice qué hacer en cada paso. Los ejemplos muestran cómo se ve ese proceso cuando los datos son reales —incluyendo las decisiones ambiguas, los modelos que no funcionan a la primera, y la iteración necesaria para llegar a un resultado defendible.

Los tres casos de este módulo representan tres tipos frecuentes de series de tiempo en la práctica: una serie con múltiples estacionalidades (ingeniería), una con tendencia y estacionalidad económica (economía), y una puramente estacional sin tendencia (ciencias ambientales). Cada una requiere decisiones de modelado diferentes.

Tres casos de uso completos con código R reproducible. Cada ejemplo sigue el flujo de los siete pasos del Módulo 6: exploración → transformación → tendencia → estacionalidad → ajuste → diagnóstico → pronóstico.

7.1 Ingeniería — Predicción de demanda eléctrica horaria

La demanda eléctrica horaria presenta múltiples estacionalidades: un patrón diario (24 horas), uno semanal (7 días) y, a veces, uno anual. Este ejemplo trabaja con datos horarios de una semana.

Características de los datos

Frecuencia: horaria ( $s = 24$ para patrón diario).
Tendencia: polinomial de bajo grado (posiblemente lineal).
Estacionalidad: doble — hora del día y día de la semana.
Tipo: generalmente aditivo (varianza aproximadamente constante).

$Y_t = \beta_0 + \beta_1 t + \sum_{i=1}^{23} \delta_i I_{i,t}^{(\text{hora})} + \sum_{j=1}^{6} \gamma_j I_{j,t}^{(\text{día})} + E_t$

Modelo horario con indicadoras de hora y día

1library(forecast)
2library(lubridate)
3
4# ── Cargar y preparar datos ───────────────────────────────
5# consumo: vector de lecturas horarias
6# fechas:  vector de POSIXct
7n   <- length(consumo)
8t   <- 1:n
9yt  <- ts(consumo, frequency = 24)  # frecuencia diaria
10
11# Variables de hora del día (ref = hora 24)
12hora_dum <- seasonaldummy(yt)  # 23 indicadoras
13
14# Variables de día de la semana (ref = domingo)
15dia_semana <- factor(weekdays(fechas), levels =
16c("lunes","martes","miércoles","jueves","viernes","sábado","domingo"))
17dia_dum <- model.matrix(~ dia_semana)[, -1]  # ref: domingo
18
19# Tendencia lineal
20poli <- data.frame(t1 = t)
21
22# ── Ajuste del modelo ─────────────────────────────────────
23datos_mod <- data.frame(yt = consumo, poli, hora_dum, dia_dum)
24modelo_elec <- lm(yt ~ ., data = datos_mod)
25cat("R² ajustado:", round(summary(modelo_elec)$adj.r.squared, 4), "\n")
26
27# ── Diagnóstico rápido ────────────────────────────────────
28e <- residuals(modelo_elec)
29Box.test(e, lag = 24, type = "Ljung-Box")  # ¿independencia?
30
31# ── Pronóstico para las próximas 24 horas ─────────────────
32t_fut      <- (n+1):(n+24)
33hora_fut   <- seasonaldummy(yt, h = 24)
34# Asumir próximo día = lunes
35dia_fut    <- matrix(c(1,0,0,0,0,0), nrow = 24, ncol = 6, byrow = TRUE)
36colnames(dia_fut) <- colnames(dia_dum)
37nuevos     <- data.frame(t1 = t_fut, hora_fut, dia_fut)
38pred_elec  <- predict(modelo_elec, newdata = nuevos, interval = "prediction")
39
40# ── Visualización ─────────────────────────────────────────
41plot(t[1:168], consumo[1:168], type = "l", col = "#d6d3d1", lwd = 1,
42   main = "Demanda eléctrica: observado vs ajustado (1 semana)",
43   xlab = "Hora", ylab = "kWh")
44lines(t[1:168], fitted(modelo_elec)[1:168], col = "#1d4ed8", lwd = 2)

Para datos con múltiples estacionalidades considere también los modelos TBATS o mstl() de la librería forecast.

Extensión a múltiples estacionalidades

Con datos horarios que abarcan varios años, la doble estacionalidad (diaria + semanal) se puede modelar con armónicos de Fourier en lugar de indicadoras, reduciendo el número de parámetros de $23 + 6 = 29$ a potencialmente 4–6 parámetros con $k = 2$ armónicos por ciclo. La función fourier() del paquete forecast facilita esto.

7.2 Economía — PIB trimestral del sector agrícola (Colombia)

Serie del PIB trimestral del sector agropecuario de Colombia en miles de millones de pesos constantes. Presenta tendencia creciente con variabilidad que aumenta con el nivel → candidato a modelo log-cuadrático.

Lección de este ejemplo

El PIB agrícola de Colombia ilustra una característica frecuente en series económicas colombianas: crecimiento con estacionalidad trimestral marcada (primer trimestre suele ser el más bajo por la temporada de siembra; cuarto trimestre el más alto por cosechas y cierre de año). Los coeficientes $\exp(\hat{\delta}_i)$ revelan directamente el porcentaje del PIB de cada trimestre relativo al cuarto trimestre de referencia —una lectura económicamente interpretable.

Análisis completo: PIB agrícola trimestral

1library(forecast)
2
3# ── Datos (reemplazar con datos reales del DANE) ──────────
4# yt: ts trimestral, ej. desde 2000-Q1
5n <- length(yt)
6t <- 1:n
7
8# ── Paso 1: Exploración ───────────────────────────────────
9plot(yt, main = "PIB sector agropecuario Colombia",
10   ylab = "Miles de millones COP (constantes)", col = "#1d4ed8")
11plot(log(yt), main = "log(PIB agropecuario)", col = "#1d4ed8")
12ggseasonplot(yt, year.labels = TRUE)  # ¿cambia la estacionalidad?
13
14# ── Paso 2: Descomposición exploratoria ───────────────────
15descom <- decompose(log(yt), type = "additive")
16plot(descom)
17
18# ── Paso 3: Modelado de tendencia ─────────────────────────
19# Observar la tendencia extraída sin los ciclos
20plot(descom$trend, main = "Tendencia estimada")
21# Grado 2 parece adecuado para el crecimiento observado
22
23poli  <- Mipoly(t, grado = 2)
24trim  <- seasonaldummy(yt)             # Q1, Q2, Q3 (ref Q4)
25
26# ── Paso 4: Ajuste del modelo log-cuadrático con dummies ──
27mod_pib <- lm(log(as.numeric(yt)) ~ ., data = data.frame(poli, trim))
28mse_pib <- summary(mod_pib)$sigma^2
29cat("\nResumen del modelo:\n")
30print(summary(mod_pib)$coefficients)
31cat("\nR² ajustado:", round(summary(mod_pib)$adj.r.squared, 4), "\n")
32
33# ── Efectos estacionales ──────────────────────────────────
34deltas   <- coef(mod_pib)[paste0("Q", 1:3)]
35efec_rel <- exp(deltas)
36cat("\nEfecto relativo vs Q4:\n")
37for (i in 1:3) {
38cat(sprintf("  Q%d: %.4f (%.2f%% del PIB de Q4)\n",
39            i, efec_rel[i], efec_rel[i]*100))
40}
41
42# ── Paso 5: Diagnóstico ───────────────────────────────────
43e <- residuals(mod_pib)
44par(mfrow = c(2,2))
45plot(t, e, type="l", main="Residuos vs t"); abline(h=0,col="red",lty=2)
46qqnorm(e); qqline(e, col="red")
47acf(e, lag.max=20, main="ACF residuos")
48pacf(e, lag.max=20, main="PACF residuos")
49par(mfrow=c(1,1))
50
51shapiro.test(e)
52Box.test(e, lag=8, type="Ljung-Box")
53
54# ── Paso 6: Pronóstico 4 trimestres ──────────────────────
55h       <- 4
56t_fut   <- (n+1):(n+h)
57nuevos  <- data.frame(Mipoly(t_fut, 2), seasonaldummy(yt, h))
58names(nuevos) <- names(coef(mod_pib))[-1]
59pred    <- predict(mod_pib, newdata=nuevos, interval="prediction", level=0.95)
60ythat_f <- exp(pred[,"fit"])   * exp(mse_pib/2)
61li_95   <- exp(pred[,"lwr"])   * exp(mse_pib/2)
62ls_95   <- exp(pred[,"upr"])   * exp(mse_pib/2)
63
64cat("\nPronóstico PIB (miles de millones COP):\n")
65trimestres_fut <- c("Q1","Q2","Q3","Q4")
66for (l in 1:h) {
67cat(sprintf("  %s: %.1f  [%.1f, %.1f]\n",
68            trimestres_fut[l], ythat_f[l], li_95[l], ls_95[l]))
69}

El código asume que la función Mipoly() está cargada en el entorno. Ajuste el grado del polinomio según el AIC.

7.3 Ambiente — Temperatura mensual de Nottingham Castle

Serie de temperatura media mensual (°F) registrada en Nottingham Castle, 1920–1939 ( $n = 240$ , $s = 12$ ). Es un ejemplo clásico de serie puramente estacional sin tendencia: la media global es aproximadamente constante.

Modelo apropiado

Sin tendencia, el modelo se reduce a:

$Y_t = \beta_0 + \sum_{i=1}^{11} \delta_i I_{i,t} + E_t$

con 12 meses y referencia en diciembre. $\beta_0$ estima la temperatura media de diciembre; $\delta_i$ la diferencia absoluta respecto a diciembre para el mes $i$ .

▶ Modelo solo estacional — Nottingham Castle (nottem)(forecast, ggplot2)

library(forecast)
library(ggplot2)

# ── Datos (disponibles en R base) ────────────────────────
yt <- nottem        # ts mensual 1920-1939, °F
n  <- length(yt)
t  <- 1:n

# ── Exploración con ggplot2 ───────────────────────────────
df_yt <- data.frame(
Año = as.numeric(time(yt)),
Tmp = as.numeric(yt)
)
mod_lineal <- lm(Tmp ~ Año, data = df_yt)
pval_tend  <- summary(mod_lineal)$coef[2, 4]
cat(sprintf("P-valor tendencia: %.4f → %s\n",
          pval_tend,
          if (pval_tend > 0.05) "sin tendencia significativa"
          else "tendencia significativa"))

ggplot(df_yt, aes(x = Año, y = Tmp)) +
geom_line(colour = "#1d4ed8", linewidth = 1) +
geom_smooth(method = "lm", se = FALSE,
            colour = "#ef4444", linetype = "dashed",
            linewidth = 0.8) +
labs(title = "Temperatura mensual — Nottingham Castle (°F)",
     subtitle = paste0("Línea roja = tendencia lineal (p = ",
                       round(pval_tend, 3), ")"),
     x = "Año", y = "°F") +
theme_bw(base_size = 12)

# ── Modelo solo estacional ─────────────────────────────────
mes_dum   <- seasonaldummy(yt)      # 11 columnas: Ene…Nov (ref = Dic)
datos_est <- data.frame(yt = as.numeric(yt), mes_dum)
mod_est   <- lm(yt ~ ., data = datos_est)
cat("\nR² ajustado:", round(summary(mod_est)$adj.r.squared, 4), "\n")

# ── Efectos mensuales ─────────────────────────────────────
coefs  <- coef(mod_est)
beta0  <- coefs[1]
deltas <- coefs[-1]
meses  <- c("Ene","Feb","Mar","Abr","May","Jun",
          "Jul","Ago","Sep","Oct","Nov")
cat(sprintf("Media de diciembre (referencia): %.2f °F\n", beta0))
for (i in seq_along(meses))
cat(sprintf("  %s: %.2f °F  (delta = %+.2f)\n",
            meses[i], beta0 + deltas[i], deltas[i]))

# ── Gráfico de efectos con ggplot2 ────────────────────────
todos_meses <- c(meses, "Dic")
temp_meses  <- c(deltas, 0) + beta0
df_bar <- data.frame(
Mes   = factor(todos_meses, levels = todos_meses),
Temp  = temp_meses,
Sobre = temp_meses >= beta0
)

ggplot(df_bar, aes(x = Mes, y = Temp, fill = Sobre)) +
geom_col(show.legend = FALSE) +
geom_hline(yintercept = beta0,
           colour = "#1c1917", linetype = "dashed", linewidth = 1) +
scale_fill_manual(values = c("TRUE" = "#ef4444", "FALSE" = "#3b82f6")) +
scale_y_continuous(limits = c(30, 70)) +
labs(title    = "Temperatura media estimada por mes (°F)",
     subtitle = paste0("Línea = referencia diciembre (",
                       round(beta0, 1), " °F)"),
     x = NULL, y = "°F") +
theme_bw(base_size = 12)

# ── Diagnóstico ───────────────────────────────────────────
e_est <- residuals(mod_est)
cat("\n--- Shapiro-Wilk ---\n"); print(shapiro.test(e_est))
cat("--- Ljung-Box ---\n")
print(Box.test(e_est, lag = 24, type = "Ljung-Box"))

# ── Pronóstico (repetir patrón estacional) ────────────────
h        <- 12
mes_fut  <- as.data.frame(seasonaldummy(yt, h = h))
pred_not <- predict(mod_est,
                  newdata  = mes_fut,
                  interval = "prediction", level = 0.95)

cat("\nProyección 1 año adelante:\n")
print(data.frame(
Mes  = todos_meses,
Pron = round(pred_not[, "fit"], 1),
LI95 = round(pred_not[, "lwr"], 1),
LS95 = round(pred_not[, "upr"], 1)
))

Ctrl+Enter

La serie nottem está disponible en R base. Es un banco de pruebas clásico para modelos estacionales sin tendencia.

Resultado clave

El mes más cálido es julio (estimado ≈ 63.5°F) y el más frío es enero (≈ 39.4°F). La diferencia de ~24°F entre julio y diciembre corresponde a $\hat{\delta}_{\text{Jul}} \approx 21$ . El modelo explica más del 97% de la variabilidad total de la serie.

7.4 Resumen comparativo de los tres ejemplos

| Aspecto | Demanda eléctrica | PIB agrícola | Temp. Nottingham | |---|---|---|---| | Frecuencia | Horaria (s=24) | Trimestral (s=4) | Mensual (s=12) | | Tendencia | Lineal leve | Cuadrática creciente | Sin tendencia | | Tipo de modelo | Aditivo | Log-cuadrático | Aditivo (solo S_t) | | Estacionalidad | Doble (hora + día) | Simple (trimestral) | Simple (mensual) | | Predictores | t, I_hora, I_día | t, t², I_Q1, I_Q2, I_Q3 | I_Ene,…,I_Nov | | Transformación | Ninguna | log(Y_t) | Ninguna |