Metodología de forecast-arena¶

Este documento explica los fundamentos estadísticos de la librería: qué hace cada pieza, por qué lo hace así, y dónde leer más. Las referencias son a Hyndman & Athanasopoulos, Forecasting: Principles and Practice, 3.ª ed. (fpp3), el libro que inspiró este proyecto.

1. Intervalos de predicción — porque un número solo es media respuesta¶

La analogía. Si un meteorólogo te dice "mañana lloverán 20 mm", te está ocultando lo importante: ¿está seguro entre 18 y 22, o entre 5 y 60? La decisión que tomas (¿llevo paraguas o cancelo el evento?) depende más del rango que del número. Con los pronósticos de negocio pasa igual: "ventas de agosto: 450" no permite decidir cuánto inventario comprar; "450, y con 95% de confianza entre 380 y 530" sí.

Qué hace la librería. Todo resultado incluye result.forecast_intervals: el pronóstico puntual más bandas al 80% y 95% (los niveles estándar de fpp3, configurables con levels=). Cada luchador calcula sus bandas con el método que le corresponde a su teoría:

Modelo	Método	Referencia
Naive	Fórmula cerrada: σₕ = σ·√h	fpp3 cap. 5.5
SNaive	σₕ = σ·√(k+1), k = ⌊(h−1)/m⌋	fpp3 cap. 5.5
Drift	σₕ = σ·√(h·(1+h/n))	fpp3 cap. 5.5
SARIMA / Harmonic	Intervalos analíticos del modelo estimado	fpp3 cap. 9.5
ETS	Distribución predictiva del espacio de estados	fpp3 cap. 8.7
XGB (híbrido)	Simulación bootstrap (ver abajo)	fpp3 cap. 5.5

Por qué el ancho crece con el horizonte. Pronosticar a 1 mes es adivinar un paso en la niebla; a 12 meses, doce pasos, y cada paso hereda la incertidumbre del anterior. Las fórmulas de arriba capturan exactamente esa acumulación: por eso las bandas se abren como un cono.

El caso especial del modelo de machine learning. XGBoost no trae fórmula de intervalos, y hay una trampa conocida: sus residuales in-sample son casi cero (el boosting memoriza el conjunto de entrenamiento), así que usarlos daría bandas de ancho cero — una mentira de precisión. La librería lo resuelve en dos partes:

Residuales honestos: se reserva internamente el último 25% del entrenamiento como holdout; los errores a un paso sobre ese tramo son la medida real de cuánto se equivoca el modelo con datos que no vio.
Bootstrap de trayectorias: se simulan 500 futuros posibles; en cada paso, la predicción recursiva se perturba con un residual honesto muestreado al azar, y ese valor perturbado se realimenta al siguiente paso. Los percentiles de las 500 trayectorias son las bandas. Así la incertidumbre se propaga por la recursión igual que lo haría en la realidad.

Intervalos y transformaciones. Cuando la Arena aplicó Box-Cox o logaritmo, las bandas se destransforman directamente: como la inversa de Box-Cox es monótona, el cuantil 97.5% en escala transformada sigue siendo el cuantil 97.5% en escala original — los intervalos son exactos. Dos consecuencias visibles y correctas: (a) las bandas quedan asimétricas (la cola superior es más larga, natural en series que crecen multiplicativamente), y (b) el punto central pasa a ser la mediana del pronóstico, no la media (fpp3 cap. 5.6). Para decisiones de inventario o presupuesto, la mediana suele ser incluso lo que quieres.

Una advertencia honesta. Cuando usas exógenas pronosticadas en cascada (forecast_X=True), los intervalos se calculan condicionales a ese pronóstico de las exógenas — es decir, no incluyen la incertidumbre extra de haber tenido que adivinar las X. Las bandas reales son algo más anchas de lo reportado. Propagar esa incertidumbre (simulando también trayectorias de las exógenas) está en el roadmap; mientras tanto, result.exog_report te dice qué tan fiable fue cada exógena para que calibres cuánto ensanchar mentalmente.

2. Diagnóstico de residuales — el examen médico del campeón¶

La analogía. Que un corredor gane la carrera no significa que esté sano. El leaderboard te dice quién ganó; el diagnóstico de residuales te dice si el ganador todavía puede dar más. Son dos preguntas distintas y fpp3 (cap. 5.4) insiste en hacer las dos.

El principio. Un residual es lo que el modelo no explicó: residual = observado − ajustado. Si el modelo extrajo toda la estructura de la serie, lo que queda debe ser ruido blanco: sin patrón, sin autocorrelación, con media cero. Si los residuales aún tienen estructura (por ejemplo, se parecen a su valor de hace 12 meses), el modelo dejó información en la mesa y el pronóstico es mejorable — quizá faltó un término estacional, una exógena, o más lags.

Qué hace la librería. Tras coronar al campeón, corre automáticamente:

Test de Ljung-Box sobre sus residuales, con el lag que recomienda el libro: min(10, n/5) para datos no estacionales, min(2m, n/5) para estacionales. La hipótesis nula es "no hay autocorrelación": un p-valor alto (> 0.05) es buena noticia — no se detecta estructura remanente.
Chequeo de sesgo: si la media de los residuales es estadísticamente distinta de cero (|media| > 2·EE), el modelo pronostica sistemáticamente de más o de menos.

El resultado vive en result.residual_diagnostics, con un veredicto legible:

"Residuales consistentes con ruido blanco (Ljung-Box p=0.886 en lag 6) y sin sesgo. El modelo extrajo la estructura disponible."

o bien:

"Atención — autocorrelación remanente (Ljung-Box p=0.011 en lag 24): el modelo dejó información en la mesa."

La letra chica que el libro subraya. Residuales limpios no garantizan buenos pronósticos (eso lo mide el backtesting, que la Arena ya hace); pero residuales sucios sí garantizan que el modelo es mejorable. Por eso el diagnóstico complementa al leaderboard, no lo sustituye. Para el modelo híbrido de ML, el diagnóstico usa los residuales honestos del holdout — los in-sample dirían "todo perfecto" por el sobreajuste del boosting, que es exactamente el tipo de autoengaño que este examen existe para evitar.

3. Estacionalidad múltiple — cuando la serie baila dos ritmos a la vez¶

La analogía. Las ventas de una tienda tienen un ritmo semanal (los sábados venden más que los martes) y encima un ritmo anual (diciembre vende más que febrero). Son dos melodías sonando al mismo tiempo. SARIMA y ETS son instrumentos de una sola melodía: les puedes pedir que atiendan el ciclo semanal o el anual, pero no ambos (fpp3 cap. 12.1).

La solución del libro: regresión armónica dinámica (fpp3 cap. 10.5 y 12.1). En lugar de pedirle al modelo que "aprenda" cada estacionalidad, se la damos dibujada: pares de ondas seno/coseno (términos de Fourier) con el periodo de cada ciclo entran como regresores. Cualquier patrón estacional suave se puede reconstruir sumando unas pocas de estas ondas — es la misma idea que descomponer un acorde en sus notas. Lo que las ondas no expliquen (la dinámica de corto plazo) lo modela un ARIMA no estacional sobre los errores. Ventaja adicional: las ondas son funciones del tiempo, así que sus valores futuros se conocen con exactitud — son la exógena perfecta.

Cómo se usa. Pasa una lista de periodos y la Arena hace el resto:

# Datos diarios: ciclo semanal (7) y mensual (30) simultáneos
arena = Arena(season_length=[7, 30])
result = arena.compete(y, test_size=14, h=14)

Al detectar la lista, la Arena suma automáticamente al luchador Harmonic (que usa todos los periodos) mientras los modelos clásicos usan el primer periodo de la lista — pon primero el dominante. Y como siempre, todos compiten: si tu serie en realidad solo tenía un ritmo, un modelo clásico ganará el leaderboard y te habrás enterado sin creerle a nadie.

El parámetro K (armónicos por periodo, default 3) controla la flexibilidad: K pequeño = estacionalidad suave, K grande = más ondulada, con tope en periodo/2. Si el diagnóstico de residuales del Harmonic reporta estructura remanente en un lag estacional, súbele a K.

Cómo se conectan las tres piezas¶

El flujo completo de la Arena quedó así:

Detección: periodicidad (una o varias) y transformación, reportadas.
Pelea: backtesting rolling-origin honesto — con exógenas en cascada si las hay — y leaderboard en escala original.
Coronación: el campeón reentrena con toda la serie.
Pronóstico completo: puntual + intervalos 80/95%.
Examen médico: Ljung-Box y sesgo sobre los residuales del campeón, con veredicto en español.

La filosofía de fondo no cambió desde la v0.1: ningún modelo se cree por default — todos pelean, los benchmarks naive siempre están en el ring, y cada afirmación de la librería (la transformación elegida, el periodo detectado, la calidad del campeón) viene acompañada de la evidencia que la sustenta.

Apéndice: los luchadores especialistas (v0.5)¶

Theta. En la competencia M3 (3,003 series reales) el método Theta venció a modelos mucho más sofisticados. Hyndman & Billah (2003) demostraron que equivale a un suavizamiento exponencial simple con drift — la mitad de su mérito es recordarnos que lo simple gana seguido. Por eso está en el lineup default: es un benchmark exigente que los modelos complejos deben justificar vencer.

Croston/SBA. Cuando el 30% o más de tus observaciones son cero (demanda intermitente: refacciones, SKUs de baja rotación), los modelos clásicos se confunden — promedian los ceros con las ventas y pronostican una demanda fantasma constante que no existe. Croston separa las dos preguntas: ¿de qué tamaño es la demanda cuando ocurre? y ¿cada cuánto ocurre? Suaviza cada una por separado y pronostica su cociente. La variante SBA (default) corrige el sesgo al alza conocido del Croston original con el factor (1 − α/2). Interpretación importante: su pronóstico plano es una tasa media de demanda por periodo, no el valor que verás en un periodo concreto.

STLF (STL + ETS). Divide y vencerás: STL separa la estacionalidad (que suele ser lo más estable de una serie), ETS pronostica lo que queda (tendencia + ruido, donde vive la dificultad), y la estacionalidad futura se repite del último ciclo. Sus intervalos heredan los del ETS desplazados por la estacionalidad — un supuesto razonable cuando el patrón estacional es estable, que es exactamente el escenario donde conviene usarlo.

Sobre el ensemble ponderado. Con ensemble_weights="inverse_error", el peso de cada modelo es proporcional al inverso de la desviación estándar de sus residuales honestos (para el modelo de ML, los del holdout interno — los in-sample mentirían). Es una heurística deliberadamente simple: los pesos óptimos exactos requerirían la matriz de covarianzas de los errores entre modelos, que con pocas ventanas de backtesting se estima peor que no estimarla (Smith & Wallis, 2009, el "forecast combination puzzle": la media simple es difícil de vencer). Por eso el default sigue siendo uniforme y el leaderboard decide si el ponderado aporta.

Auditoría de intervalos (v0.6) — porque prometer 95% obliga a cumplir 95%¶

La analogía. Un pronosticador del clima que dice "95% de probabilidad de lluvia" 100 veces debe acertar unas 95. Si solo acierta 70, sus probabilidades están mal calibradas aunque suene confiado. Con las bandas de predicción es igual: decir "intervalo del 95%" es una promesa auditable.

Qué se mide (fpp3 cap. 5.9), en cada ventana del backtesting y en escala original:

Cobertura empírica: fracción de observaciones reales dentro de la banda. Debe rondar el nivel nominal. Sub-cobertura = bandas que mienten angosto (el peligro real); sobre-cobertura = bandas inútilmente anchas.
Winkler score: ancho de la banda + penalización 2/α por cada observación que escapa. Es la métrica anti-trampa: bandas infinitas cubren todo pero pagan su ancho; bandas de ancho cero son angostísimas pero pagan cada fallo. Menor = mejor.

Lo que la auditoría destapó al estrenarse (y por qué existe): las bandas bootstrap del modelo híbrido sub-cubrían (67% real en la banda "del 95%") porque les faltaba la incertidumbre de la extrapolación de tendencia. Se corrigió agregando ese término (con la forma del drift: crece como √(h·(1+h/n))) y la cobertura subió a 83%. La banda del SNaive en series con tendencia cubre casi 0% — no es un bug: sus intervalos son honestos al modelo, y el modelo está sesgado; Winkler lo castiga con números enormes y el leaderboard ya lo tenía al fondo. La auditoría convierte estas historias en una tabla que ves en cada corrida.

Pronóstico jerárquico (v0.6) — que las cuentas cuadren¶

El problema. Pronosticas cada tienda, cada región y el total nacional por separado. Cada pronóstico individual puede ser bueno, pero la suma de las tiendas NO da el total pronosticado — y en la junta de planeación, dos áreas trabajando con niveles distintos traen números incompatibles.

La solución (fpp3 cap. 11): definir la jerarquía como árbol y reconciliar. La matriz de suma S codifica quién suma a quién; reconciliar es proyectar los pronósticos base al espacio donde las sumas cuadran:

bottom_up: solo se cree a las hojas; los agregados son su suma. Bueno cuando las hojas tienen señal propia clara.
top_down: solo se cree al total; se reparte con las proporciones históricas promedio. Robusto arriba, ciego a la dinámica de cada hoja.
ols (default): ỹ = S(SᵀS)⁻¹Sᵀŷ — usa la información de TODOS los niveles a la vez y encuentra la corrección mínima que hace coherente al conjunto. Es el caso identidad-ponderado del método MinT (cap. 11.4) y en la práctica suele mejorar la precisión además de cuadrar las cuentas: los niveles agregados (menos ruidosos) corrigen a las hojas y viceversa.

Cada nodo corre su propia Arena — con su campeón, su transformación y su diagnóstico — y hr.summary te muestra el panorama completo.

Efectos de calendario (v0.6) — el ruido con fecha de caducidad conocida¶

Dos efectos que fpp3 (cap. 10.2) destaca porque contaminan la estacionalidad si no se modelan, y que son gratis de predecir (son funciones puras del calendario):

Días hábiles (trading_days): un marzo con 23 días hábiles vende más que un febrero con 20 aunque la demanda diaria sea idéntica. Sin esta variable, el modelo confunde "efecto marzo" con "efecto tres días extra".
Semana Santa (is_easter): la fiesta móvil por excelencia — cae en marzo o abril según el año. Ningún patrón estacional fijo puede capturarla: un año infla marzo, el siguiente infla abril, y el modelo ve ruido donde hay una causa perfectamente predecible.

Ambas van directo en X_future (se calculan para cualquier fecha futura): son el ejemplo canónico de exógena conocida.

Refinamientos de honestidad estadística (v0.6)¶

Propagación de la incertidumbre de la cascada. Cuando las exógenas se pronostican (forecast_X=True), las bandas del modelo son condicionales a ese pronóstico. Ahora se corrige por método delta: cada exógena se perturba ±σ (el error medido en su holdout), se mide cuánto mueve al pronóstico de y, y esa varianza extra ensancha las bandas: half² → half² + z²·Σδ². Es una aproximación de primer orden — exacta para efectos lineales (SARIMAX), aproximada para el híbrido.

Ajuste de sesgo Box-Cox (bias_adjust=True, cap. 5.6). La destransformación directa entrega la mediana del pronóstico. Para decisiones "de percentil" (¿cuánto inventario para cubrir la demanda típica?) la mediana es correcta. Pero cuando los pronósticos deben sumar (presupuestos, totales por región), se necesita la media, y media ≠ mediana en distribuciones asimétricas. El ajuste aplica E[y] ≈ w⁻¹(μ)·[1 + σ_h²(1−λ)/(2(λμ+1)²)] — siempre ≥ la mediana en series positivas, y más grande cuanto más incierto el horizonte.

K por AICc en Harmonic. El número de armónicos se elige minimizando AICc sobre una regresión rápida (Fourier + tendencia), el criterio del cap. 10.5. Ajustar el ARIMA completo por cada K sería carísimo; la aproximación OLS ordena a los candidatos casi igual y es instantánea.

Detección de quiebres estructurales (v0.7) — cuando la serie cambia de vida¶

La analogía. Pídele a alguien que pronostique tu gasto mensual usando también los años en que vivías con tus papás: esos datos son reales, pero describen una vida que ya no es la tuya, y contaminan la respuesta. Las series de negocio viven lo mismo — un cambio de precios, la apertura de un competidor, una remodelación — y un modelo que no sabe del quiebre pronostica el mundo que ya no existe.

El algoritmo: PELT (Pruned Exact Linear Time; Killick, Fette & Eckley, 2012). Encuentra la partición óptima de la serie en segmentos internamente coherentes, balanceando dos fuerzas: cada segmento debe explicarse bien por sí solo, pero cada quiebre nuevo paga una penalización. Es exacto (no heurístico), rápido gracias a su poda, y es el mismo motor de los paquetes especializados ruptures (Python) y changepoint (R). Aquí está implementado internamente, sin dependencias nuevas.

Decisiones de diseño que importan (y las lecciones de calibrarlo):

Costo piecewise-lineal, no de medias. Un detector de saltos de media sobre una serie con tendencia la convierte en "escalera": reporta un quiebre falso cada pocos meses porque una rampa se explica como muchos escalones. Lo vimos en la calibración (7-8 quiebres fantasma en series perfectamente limpias). La solución: el detector ajusta rectas por segmento — una tendencia limpia es UNA recta, cero quiebres — y después clasifica cada quiebre: si al empatar las rectas de ambos lados hay discontinuidad grande (>3σ), la serie brincó → nivel; si las rectas casi se tocan pero la pendiente cambió → tendencia.
Sobre la serie desestacionalizada (STL robusto), para que el ciclo normal no dispare alarmas.
Ruido estimado robustamente (MAD de las diferencias): la estimación de "cuánto tiembla normalmente la serie" debe ser inmune a los propios quiebres y outliers que busca.
Penalización tipo BIC, deliberadamente conservadora: preferimos no ver un quiebre sutil (2σ) a inundarte de fantasmas. La calibración final: 0 falsos positivos en 60 series limpias, 100% de detección en quiebres francos con localización casi exacta. penalty_scale ajusta la sensibilidad si tu caso pide otra cosa.

Las tres reacciones. Detectar es la mitad; qué hacer es la otra:

report: la opción epistémicamente humilde — la Arena no puede saber si el quiebre de tu serie fue un cambio permanente o un evento que se revertirá; te lo muestra y decides.
trim: amputar la historia pre-quiebre. Correcto cuando el pasado es activamente engañoso; costoso porque tira datos (por eso la guarda de historia mínima que rehúsa dejarte con media temporada).
dummy: la vía del cap. 10.4 de fpp3 — el quiebre entra como variable de intervención (escalón 0→1). El modelo conserva toda la muestra pero con la nota "aquí cambió el mundo": estima el tamaño del salto como un coeficiente más y pronostica desde el régimen nuevo. Suele dominar a las otras dos porque no desperdicia información.

La letra chica. Ningún detector distingue un quiebre estructural de un outlier sostenido o del inicio de una curva suave — matemáticamente pueden ser idénticos en muestra corta. Por eso el default es reportar, la gráfica marca los quiebres para tu inspección visual, y la decisión de fondo (¿esto fue permanente?) sigue siendo del analista que conoce el negocio.