Módulo 9·Medias móviles, Henderson, SES, Holt-Winters y Descomp+LOESS

Modelos y Métodos Locales de Suavizamiento

Filtros lineales (medias móviles unilaterales/bilaterales, Henderson), suavizamiento exponencial simple (SES), Holt-Winters estacional (SEHW) y descomposición combinada con LOESS (DLL/DLC). Criterios AICC, GCV, validación cruzada y medidas de precisión de pronósticos.

Medias móviles (unilateral y bilateral)Filtros de HendersonSuavizamiento exponencial simple (SES)Holt-Winters estacional (SEHW)Descomposición + LOESS (DLL/DLC)Criterios AICC y GCVRMSE, MAPE, Interval Score

Anterior: Pruebas de Hipótesis y Gráficos Diagnósticos Siguiente: Procesos Estocásticos Estacionarios y Funciones de Autocorrelación

Del modelo global al modelo local

Los módulos 2–4 presentaron modelos globales: una única ecuación de regresión (polinomial, log-polinomial o exponencial) ajustada por mínimos cuadrados sobre toda la serie. El módulo 5 introdujo LOESS y STL como alternativas no paramétricas. Este módulo completa el panorama de métodos locales con una familia de técnicas que estiman la tendencia y la estacionalidad usando ventanas móviles o pesos exponenciales decrecientes, sin imponer una forma funcional global.

La idea común es: el valor ajustado en el tiempo $t$ depende principalmente de las observaciones cercanas a $t$ , no de toda la serie. Esto permite capturar tendencias y patrones estacionales que cambian a lo largo del tiempo.

Medias móviles, filtros de Henderson, suavizamiento exponencial (SES y Holt-Winters estacional) y descomposición combinada con LOESS. Cada método se presenta con su fundamento teórico, código R reproducible y criterios de comparación (AIC, BIC, RMSE, MAPE, amplitud y cobertura de intervalos).

9.1 Medias móviles

La media móvil de orden $m+1$ es el promedio aritmético de $m+1$ observaciones consecutivas. Se trata de un filtro lineal: el valor ajustado en $t$ es una combinación lineal de observaciones vecinas con pesos $w_i$ .

$\hat{Y}_t = \sum_{i=-m}^{m} w_i \, Y_{t+i}$

Los pesos $w_i$ satisfacen $\sum w_i = 1$ . Según la simetría de los pesos, el filtro es:

Unilateral (hacia atrás)

Usa $Y_{t}, Y_{t-1}, \ldots, Y_{t-m}$ . Solo información pasada y presente. Útil para pronóstico en tiempo real.

\hat{Y}_t = \frac{1}{m+1}\sum_{i=0}^{m} Y_{t-i}

Bilateral (centrada)

Usa $Y_{t-m}, \ldots, Y_{t}, \ldots, Y_{t+m}$ . Promedia información pasada y futura. Suaviza mejor pero no pronostica.

\hat{Y}_t = \frac{1}{2m+1}\sum_{i=-m}^{m} Y_{t+i}

Estrategia circular para datos estacionales

Cuando la serie tiene período estacional $s$ , conviene elegir $m+1 = s$ (unilateral) o $2m+1 = s$ (bilateral) para que el filtro promedie exactamente un ciclo completo y elimine la estacionalidad de la tendencia estimada. Con $s = 12$ (mensual), se usan $m+1 = 13$ (unilateral) o $2m+1 = 25$ (bilateral).

▶ Medias móviles unilaterales y bilaterales — datos MSFT(forecast)

Ctrl+Enter

La función ma() del paquete forecast implementa medias móviles. centre=TRUE produce una media centrada (bilateral).

9.2 Filtros de Henderson

Los filtros de Henderson son medias móviles bilaterales con pesos no uniformes, diseñados para seguir tendencias suavemente (reducir el "zigzag" que producen las medias móviles simples). Los pesos se derivan minimizando la suma de los cuadrados de las terceras diferencias de la serie suavizada, sujeto a restricciones de insesgamiento para polinomios de grado ≤ 3.

Pesos de Henderson

Para un filtro de $2m+1$ términos, los pesos $w_i$ ( $i = -m, \ldots, 0, \ldots, m$ ) son:

$w_i = \left[ 315 \left((m+1)^2 - i^2\right)^2 - 21 \left((m+1)^2 - i^2\right) \left(3m(m+2) - 7\right) \right] / K$

donde $K$ es una constante de normalización tal que $\sum w_i = 1$ . Los filtros más usados son de 13 términos ( $m = 6$ ) y 23 términos ( $m = 11$ ).

Los pesos de Henderson son simétricos y suman 1. Su forma es de campana: pesos máximos en el centro ( $i=0$ ) y decrecientes hacia los extremos, con valores ligeramente negativos en los bordes para mejorar el seguimiento de curvas.

Filtro de Henderson 13 y 23 términos — datos MSFT

1w13 <- c(-0.019, -0.028, 0.0, 0.066, 0.147, 0.214, 0.240,
2        0.214, 0.147, 0.066, 0.0, -0.028, -0.019)
3sum(w13)
4w23 <- c(-0.004, -0.011, -0.016, -0.015, -0.005, 0.013, 0.039,
5        0.068, 0.097, 0.122, 0.138, 0.148, 0.138, 0.122, 0.097,
6        0.068, 0.039, 0.013, -0.005, -0.015, -0.016, -0.011, -0.004)
7sum(w23)
8fh13 <- filter(yt, filter = w13, sides = 2)
9fh23 <- filter(yt, filter = w23, sides = 2)
10plot(yt, main = "Filtros de Henderson", col = "gray60")
11lines(fh13, col = "#059669", lwd = 2)
12lines(fh23, col = "#7c3aed", lwd = 2)
13legend("topright", legend = c("Observado", "Henderson 13", "Henderson 23"),
14     col = c("gray60", "#059669", "#7c3aed"), lwd = c(1, 2, 2))
15rmse_fh13 <- sqrt(mean((yt - fh13)^2, na.rm = TRUE))
16rmse_fh23 <- sqrt(mean((yt - fh23)^2, na.rm = TRUE))
17cat(sprintf("RMSE Henderson 13: %.2f\n", rmse_fh13))
18cat(sprintf("RMSE Henderson 23: %.2f\n", rmse_fh23))

Los pesos de Henderson son fijos (no se estiman). La función filter() de R base aplica el filtro lineal.

9.3 Suavizamiento exponencial simple (SES)

El SES es el método local más parsimonioso: el valor ajustado en $t+1$ es un promedio ponderado entre la observación actual y el valor ajustado anterior, con peso $\alpha$ para la observación nueva y $1-\alpha$ para el ajuste previo.

$\hat{Y}_{t+1} = \alpha \, Y_t + (1 - \alpha) \, \hat{Y}_t$

Al iterar hacia atrás, $\hat{Y}_{t+1}$ es un promedio ponderado de todas las observaciones pasadas con pesos que decaen exponencialmente:

$\hat{Y}_{t+1} = \alpha \sum_{j=0}^{t-1} (1-\alpha)^j \, Y_{t-j} + (1-\alpha)^t \, \hat{Y}_1$

El parámetro $\alpha \in (0,1)$ controla la rapidez de adaptación: $\alpha \to 1$ da casi todo el peso a la observación más reciente (adaptación rápida); $\alpha \to 0$ promedia con pesos casi uniformes (suavizamiento fuerte). El valor óptimo se estima minimizando la suma de cuadrados de los errores de un paso adelante.

Cuándo usar SES

Series sin tendencia ni estacionalidad aparentes.
Series de alta frecuencia (diaria, horaria) donde los modelos paramétricos son costosos.
Como baseline para comparar contra modelos más complejos.

▶ SES óptimo — datos MSFT(forecast)

Ctrl+Enter

HoltWinters() con beta=FALSE y gamma=FALSE desactiva los componentes de tendencia y estacionalidad, dejando SES puro.

9.4 Suavizamiento estacional Holt-Winters (SEHW)

El método Holt-Winters extiende el SES incorporando una componente de tendencia y una componente estacional. La versión estacional (SEHW) es la más usada para series con patrón estacional periódico.

Existen dos variantes según la naturaleza de la estacionalidad:

Aditiva

Y_t = L_t + T_t + S_t + \varepsilon_t

La amplitud estacional es constante a lo largo del tiempo.

Multiplicativa

Y_t = (L_t + T_t) \cdot S_t + \varepsilon_t

La amplitud estacional crece con el nivel de la serie.

Las ecuaciones de recurrencia para la variante multiplicativa son:

$\begin{aligned} L_t &= \alpha \, \frac{Y_t}{S_{t-s}} + (1-\alpha)(L_{t-1} + T_{t-1}) \\\\ T_t &= \beta \, (L_t - L_{t-1}) + (1-\beta) \, T_{t-1} \\\\ S_t &= \gamma \, \frac{Y_t}{L_t} + (1-\gamma) \, S_{t-s} \end{aligned}$

donde $L_t$ es el nivel, $T_t$ es la tendencia, $S_t$ es el factor estacional, $s$ es el período estacional, y $\alpha, \beta, \gamma \in (0,1)$ son los parámetros de suavizamiento (estimados por máxima verosimilitud o minimizando SSE).

Número de parámetros del SEHW

El SEHW tiene $(s-1) + 2$ parámetros efectivos: $s-1$ factores estacionales (uno por período, el último se fija por normalización) más $\alpha$ y $\beta$ . El valor $\gamma$ se cuenta dentro de los $s-1$ factores. Esta aproximación se usa para calcular AIC y BIC.

SEHW multiplicativo — Ventas Cía. X (mensual, s=12)

1library(forecast)
2ventas <- c(
3154, 96, 73, 49, 36, 59, 95, 169, 210, 278, 298, 245,
4200, 118, 90, 79, 78, 91, 167, 169, 289, 347, 375, 203,
5223, 104, 107, 85, 75, 99, 135, 211, 335, 460, 488, 326,
6346, 261, 224, 141, 148, 145, 223, 272, 445, 560, 612, 467,
7518, 404, 300, 210, 196, 186, 247, 343, 464, 680, 711, 610,
8613, 392, 273, 322, 189, 257, 324, 404, 677, 858, 895, 664,
9628, 308, 324, 248, 272
10)
11ventas <- ts(ventas, frequency = 12, start = c(1965, 1))
12m <- 6; n <- length(ventas) - m
13yt <- ts(ventas[1:n], frequency = 12, start = c(1965, 1))
14ytf <- ts(ventas[(n+1):length(ventas)], frequency = 12, start = c(1970, 12))
15sehw <- HoltWinters(yt, seasonal = "multiplicative")
16sehw
17cat("Alpha:", round(sehw$alpha, 4), "Beta:", round(sehw$beta, 4), "Gamma:", round(sehw$gamma, 4), "\n")
18npar <- 12 - 1 + 2; sse <- sehw$SSE
19aic_sehw <- n * log(sse / n) + 2 * npar
20bic_sehw <- n * log(sse / n) + npar * log(n)
21cat(sprintf("AIC ≈ %.2f  BIC ≈ %.2f\n", aic_sehw, bic_sehw))
22pron_sehw <- predict(sehw, n.ahead = m, prediction.interval = TRUE)
23ytpron <- ts(pron_sehw[, 1], frequency = 12, start = c(1970, 12))
24accuracy(ytpron, ytf)
25plot(ytf, type = "b", pch = 19, lwd = 2, ylim = range(ytf, pron_sehw),
26   main = "Pronóstico SEHW vs real (Dic 1970 – May 1971)")
27lines(ytpron, col = 2, pch = 2, lty = 2, type = "b", lwd = 2)
28lines(pron_sehw[, 2], col = 2, lty = 3)
29lines(pron_sehw[, 3], col = 2, lty = 3)
30legend("topright", legend = c("Real", "Pronóstico SEHW", "IP 95%"),
31     col = c(1, 2, 2), lty = c(1, 2, 3), pch = c(19, 2, NA), lwd = 2)

HoltWinters() con seasonal='multiplicative' estima alpha, beta y gamma automáticamente.

9.5 Descomposición combinada con LOESS (DLL/DLC)

Este método, desarrollado por Nelfi González, combina un filtro de descomposición clásico (como decompose()) para extraer el componente estacional $S_t$ , con un ajuste LOESS para la tendencia $T_t$ de la serie desestacionalizada. La ventaja es que la tendencia se modela localmente (no se impone un polinomio global), mientras que la estacionalidad se obtiene con un filtro robusto.

Esquema del método Descomp+LOESS

Paso 1: Aplicar filtro de descomposición para obtener $\hat{S}_t$ :

$Y_t = T_t \times S_t \times E_t \quad \text{(multiplicativo)}$

$Y_t = T_t + S_t + E_t \quad \text{(aditivo)}$

Paso 2: Desestacionalizar: $Y_t^{des} = Y_t / \hat{S}_t$ (multiplicativo) o $Y_t^{des} = Y_t - \hat{S}_t$ (aditivo).

Paso 3: Ajustar LOESS de grado $g$ a $Y_t^{des}$ para obtener $\hat{T}_t$ .

Paso 4: Valores ajustados: $\hat{Y}_t = \hat{T}_t \times \hat{S}_t$ (multiplicativo).

El parámetro de suavizamiento $\alpha$ (span) del LOESS se selecciona automáticamente minimizando un criterio: AICC (criterio de información de Akaike corregido) o GCV (validación cruzada generalizada).

Se denotan:

DLL — Descomposición + LOESS Lineal (grado $g = 1$ )
DLC — Descomposición + LOESS Cuadrático (grado $g = 2$ )

Cada combinación se evalúa con ambos criterios (AICC y GCV), resultando en cuatro modelos: DLL(AICC), DLL(GCV), DLC(AICC), DLC(GCV).

Filtro de descomposición clásico — Ventas Cía. X

1descom <- decompose(ventas, type = "multiplicative")
2St <- descom$seasonal; Tt <- descom$trend; Et <- descom$random
3ytd <- ventas / St
4aux <- split(as.numeric(ytd), cycle(ytd))
5medias_mes <- sapply(aux, mean, na.rm = TRUE)
6deltas_est <- medias_mes / mean(medias_mes)
7cat("Efectos estacionales (multiplicativos):\n")
8print(round(deltas_est, 4))
9descom_log <- decompose(log(ventas), type = "additive")
10St_log <- descom_log$seasonal; Tt_log <- descom_log$trend; Et_log <- descom_log$random
11aux_log <- split(as.numeric(log(ventas) - Tt_log), cycle(ventas))
12medias_log <- sapply(aux_log, mean, na.rm = TRUE)
13deltas_log <- medias_log - mean(medias_log)
14cat("\nEfectos estacionales (log-escala aditiva):\n")
15print(round(deltas_log, 4))

decompose() usa un filtro de medias móviles centrado para extraer T_t y promedios por período para S_t.

Descomp+LOESS (DLL/DLC) — Ventas Cía. X

1library(forecast); library(fANCOVA)
2source("https://raw.githubusercontent.com/NelfiGonzalez/Funciones-de-Usuario-Estadistica-III/main/Funcion-Descomp.Loessv02.R")
3source("https://raw.githubusercontent.com/NelfiGonzalez/Funciones-de-Usuario-Estadistica-III/main/Funciones-Criterios.Informacion-Calidad.Intervalos.R")
4ajuste_DLL_aicc <- Descomp.Loessv02(serie.ajuste = yt, tipo.descomp = "multiplicative", grado = 1, criterio = "aicc", h = m, level = 0.95)
5ajuste_DLL_gcv <- Descomp.Loessv02(serie.ajuste = yt, tipo.descomp = "multiplicative", grado = 1, criterio = "gcv", h = m, level = 0.95)
6ajuste_DLC_aicc <- Descomp.Loessv02(serie.ajuste = yt, tipo.descomp = "multiplicative", grado = 2, criterio = "aicc", h = m, level = 0.95)
7ajuste_DLC_gcv <- Descomp.Loessv02(serie.ajuste = yt, tipo.descomp = "multiplicative", grado = 2, criterio = "gcv", h = m, level = 0.95)
8criterios <- rbind(
9exp.crit.inf.resid(residuals(ajuste_DLL_aicc), ajuste_DLL_aicc$p),
10exp.crit.inf.resid(residuals(ajuste_DLL_gcv), ajuste_DLL_gcv$p),
11exp.crit.inf.resid(residuals(ajuste_DLC_aicc), ajuste_DLC_aicc$p),
12exp.crit.inf.resid(residuals(ajuste_DLC_gcv), ajuste_DLC_gcv$p))
13rownames(criterios) <- c("DLL(AICC)", "DLL(GCV)", "DLC(AICC)", "DLC(GCV)")
14criterios
15pron_DLL_aicc <- ajuste_DLL_aicc$forecast[, 1]
16pron_DLL_gcv  <- ajuste_DLL_gcv$forecast[, 1]
17pron_DLC_aicc <- ajuste_DLC_aicc$forecast[, 1]
18pron_DLC_gcv  <- ajuste_DLC_gcv$forecast[, 1]
19comp <- rbind(accuracy(pron_DLL_aicc, ytf), accuracy(pron_DLL_gcv, ytf),
20            accuracy(pron_DLC_aicc, ytf), accuracy(pron_DLC_gcv, ytf))[, c(2, 3, 5)]
21rownames(comp) <- c("DLL(AICC)", "DLL(GCV)", "DLC(AICC)", "DLC(GCV)")
22comp

La función Descomp.Loessv02() combina el filtro de descomposición con LOESS. El parámetro 'grado' controla el grado del polinomio local de LOESS.

9.6 Comparación de métodos locales

| Método | Tendencia | Estacionalidad | Parámetros | Pronóstico | |---|---|---|---|---| | MM unilateral | Sí (local) | No | 1 (ventana) | Sí (naive) | | MM bilateral | Sí (local) | No | 1 (ventana) | No | | Henderson | Sí (suave) | No | 1 (ventana) | No | | SES | No | No | 1 (α) | Sí | | SEHW | Sí (lineal) | Sí (s períodos) | 3 + s−1 | Sí | | DLL/DLC | Sí (LOESS) | Sí (filtro) | 2 (grado, α) | Sí |

Criterios de selección

Sin tendencia ni estacionalidad: SES es la opción más parsimoniosa.
Con tendencia, sin estacionalidad: medias móviles o Henderson (para suavizar); Holt lineal para pronosticar.
Con estacionalidad estable: SEHW aditivo o DLL/DLC aditivo.
Con estacionalidad creciente: SEHW multiplicativo o DLL/DLC multiplicativo.
Tendencia no lineal cambiante: DLC (LOESS cuadrático) con criterio AICC.

9.7 Medidas de precisión de pronósticos

Todos los métodos locales se comparan mediante validación cruzada: se reservan los últimos $m$ períodos como muestra de prueba y se ajusta el modelo con los $n-m$ restantes. Las medidas de precisión son:

RMSE

\text{RMSE} = \sqrt{\frac{1}{m}\sum_{t=n+1}^{n+m} (Y_t - \hat{Y}_t)^2}

MAE

\text{MAE} = \frac{1}{m}\sum_{t=n+1}^{n+m} |Y_t - \hat{Y}_t|

MAPE

\text{MAPE} = \frac{100}{m}\sum_{t=n+1}^{n+m} \left|\frac{Y_t - \hat{Y}_t}{Y_t}\right|

Interval Score

S = (LSP - LIP) + \frac{2}{\alpha}(LIP - Y_t)\mathbb{1}_{Y_t < LIP} + \frac{2}{\alpha}(Y_t - LSP)\mathbb{1}_{Y_t > LSP}

Para los intervalos de predicción (IP) se reportan además la amplitud media $(LSP - LIP)$ y la cobertura empírica (proporción de observaciones que caen dentro del IP). Un IP del 95% bien calibrado debe tener cobertura ≈ 95%.

Tabla comparativa final — datos MSFT

1tabla_rmse <- rbind(RMSE_MM_unilateral_13, RMSE_MM_bilateral_25,
2RMSE_Henderson_13, RMSE_Henderson_23, RMSE_SES)
3rownames(tabla_rmse) <- c("MMSU13", "MMSB25", "FH13", "FH23", "SES")
4colnames(tabla_rmse) <- "RMSE de ajuste"
5tabla_rmse
6comp_pr <- rbind(
7accuracy(ytpron_MMSU13, ytnuevo), accuracy(ytpron_MMSB25, ytnuevo),
8accuracy(ytpron_FH13, ytnuevo), accuracy(ytpron_FH23, ytnuevo),
9accuracy(ytpron_SES, ytnuevo))[, c(2, 3, 5)]
10rownames(comp_pr) <- c("MMSU13", "MMSB25", "FH13", "FH23", "SES")
11colnames(comp_pr) <- c("RMSE", "MAE", "MAPE")
12comp_pr
13amp_cob <- rbind(ampcob_MMSU13, ampcob_MMSB25, ampcob_FH13, ampcob_FH23, ampcob_SES)
14score_ip <- rbind(scoreip_MMSU13, scoreip_MMSB25, scoreip_FH13, scoreip_FH23, scoreip_SES)
15tabla_final <- cbind(tabla_rmse, comp_pr, amp_cob, score_ip)
16rownames(tabla_final) <- c("MMSU13", "MMSB25", "FH13", "FH23", "SES")
17tabla_final

La función accuracy() del paquete forecast calcula RMSE, MAE y MAPE.

Conexión con TimeSight

El software TimeSight incorpora el flujo de análisis de este módulo en su paso de Modelado. Cuando selecciones la familia ARIMA automático, TimeSight usa auto.arima() que —al igual que los métodos locales— selecciona automáticamente los parámetros óptimos. Los módulos 10 y 11 profundizan en la teoría de los procesos ARMA que sustentan esta opción.