Series Temporales¶

Una serie temporal es una secuencia de observaciones ordenadas en el tiempo. El orden importa: a diferencia del resto de los datos que hemos trabajado, aquí las observaciones pasadas contienen información sobre las futuras. Este notebook cubre el análisis exploratorio de series temporales, la descomposición de sus componentes, y dos enfoques de pronóstico: el clásico con statsmodels (SARIMA) y el moderno con features de ML (scikit-learn).

Librerías¶

In [25]:

Copied!





import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import warnings
warnings.filterwarnings('ignore')

# Estadística y series temporales
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX

# ML
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

plt.rcParams['figure.figsize'] = (12, 4)
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False

np.random.seed(42)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import warnings
warnings.filterwarnings('ignore')

# Estadística y series temporales
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX

# ML
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

plt.rcParams['figure.figsize'] = (12, 4)
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False

np.random.seed(42)

1. Dataset: Ventas Mensuales¶

Trabajamos con un dataset de ventas mensuales de una cadena de retail. La serie tiene tendencia creciente, estacionalidad anual (peaks en diciembre/enero) y ruido.

In [26]:

Copied!





# Serie temporal sintética pero realista: 7 años de ventas mensuales
np.random.seed(42)
fechas = pd.date_range(start='2017-01-01', periods=84, freq='MS')  # 7 años

t = np.arange(84)
tendencia   = 800_000 + 8_000 * t                                   # crecimiento lineal
estacional  = 200_000 * np.sin(2 * np.pi * t / 12 - np.pi / 2)     # ciclo anual
peak_dic    = 300_000 * (np.array([1 if (i % 12) == 11 else 0       # peak diciembre
              for i in range(84)]))
ruido       = np.random.normal(0, 80_000, 84)

ventas = tendencia + estacional + peak_dic + ruido

serie = pd.Series(ventas, index=fechas, name='ventas')

print(f"Período: {fechas[0].date()} a {fechas[-1].date()}")
print(f"Observaciones: {len(serie)}")
print(f"Media: {serie.mean():,.0f} CLP")
print(f"Desv. estándar: {serie.std():,.0f} CLP")
# Serie temporal sintética pero realista: 7 años de ventas mensuales
np.random.seed(42)
fechas = pd.date_range(start='2017-01-01', periods=84, freq='MS')  # 7 años

t = np.arange(84)
tendencia   = 800_000 + 8_000 * t                                   # crecimiento lineal
estacional  = 200_000 * np.sin(2 * np.pi * t / 12 - np.pi / 2)     # ciclo anual
peak_dic    = 300_000 * (np.array([1 if (i % 12) == 11 else 0       # peak diciembre
              for i in range(84)]))
ruido       = np.random.normal(0, 80_000, 84)

ventas = tendencia + estacional + peak_dic + ruido

serie = pd.Series(ventas, index=fechas, name='ventas')

print(f"Período: {fechas[0].date()} a {fechas[-1].date()}")
print(f"Observaciones: {len(serie)}")
print(f"Media: {serie.mean():,.0f} CLP")
print(f"Desv. estándar: {serie.std():,.0f} CLP")

Período: 2017-01-01 a 2023-12-01
Observaciones: 84
Media: 1,148,609 CLP
Desv. estándar: 267,280 CLP

In [27]:

Copied!





# Visualización inicial
fig, ax = plt.subplots(figsize=(13, 4))
ax.plot(serie.index, serie/1e6, color='steelblue', lw=1.8)
ax.fill_between(serie.index, serie/1e6, alpha=0.1, color='steelblue')
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax.set_title('Ventas mensuales 2017–2023')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
plt.tight_layout()
plt.show()
# Visualización inicial
fig, ax = plt.subplots(figsize=(13, 4))
ax.plot(serie.index, serie/1e6, color='steelblue', lw=1.8)
ax.fill_between(serie.index, serie/1e6, alpha=0.1, color='steelblue')
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax.set_title('Ventas mensuales 2017–2023')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
plt.tight_layout()
plt.show()

No description has been provided for this image

2. Análisis Exploratorio de Series Temporales¶

Antes de modelar, es fundamental entender la estructura de la serie. Las cuatro preguntas clave son:

Pregunta	¿Qué buscar?
¿Hay tendencia?	¿La serie sube o baja sistemáticamente?
¿Hay estacionalidad?	¿Se repiten patrones en períodos fijos?
¿Hay outliers?	¿Hay valores extremos o eventos atípicos?
¿Es estacionaria?	¿Media y varianza son constantes en el tiempo?

In [28]:

Copied!





# Estadísticas por año
df_anual = serie.resample('Y').agg(['mean', 'std', 'min', 'max'])
df_anual.index = df_anual.index.year
df_anual.columns = ['Media', 'Std', 'Mínimo', 'Máximo']
print("=== Estadísticas anuales ===")
print((df_anual / 1e6).round(2).to_string())
# Estadísticas por año
df_anual = serie.resample('Y').agg(['mean', 'std', 'min', 'max'])
df_anual.index = df_anual.index.year
df_anual.columns = ['Media', 'Std', 'Mínimo', 'Máximo']
print("=== Estadísticas anuales ===")
print((df_anual / 1e6).round(2).to_string())

=== Estadísticas anuales ===
      Media   Std  Mínimo  Máximo
2017   0.89  0.17    0.62    1.17
2018   0.92  0.18    0.58    1.18
2019   1.05  0.18    0.75    1.37
2020   1.13  0.20    0.77    1.39
2021   1.26  0.20    0.88    1.51
2022   1.36  0.21    1.04    1.62
2023   1.44  0.16    1.08    1.63

In [29]:

Copied!





# Patrón estacional: ventas promedio por mes
df_mensual = pd.DataFrame({'ventas': serie.values, 'mes': serie.index.month,
                            'año': serie.index.year})
meses_nombre = ['Ene','Feb','Mar','Abr','May','Jun',
                'Jul','Ago','Sep','Oct','Nov','Dic']

fig, axes = plt.subplots(1, 2, figsize=(14, 4))

# Promedio por mes
ax = axes[0]
prom_mes = df_mensual.groupby('mes')['ventas'].mean()
colores_mes = ['crimson' if m == 12 else 'steelblue' for m in range(1, 13)]
ax.bar(range(1, 13), prom_mes/1e6, color=colores_mes, edgecolor='white', alpha=0.85)
ax.set_xticks(range(1, 13))
ax.set_xticklabels(meses_nombre)
ax.set_title('Ventas promedio por mes (2017–2023)')
ax.set_ylabel('Ventas (millones CLP)')

# Boxplot por mes
ax = axes[1]
data_box = [df_mensual[df_mensual['mes'] == m]['ventas'].values / 1e6
            for m in range(1, 13)]
bp = ax.boxplot(data_box, patch_artist=True, notch=False)
for patch, color in zip(bp['boxes'], colores_mes):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)
ax.set_xticklabels(meses_nombre)
ax.set_title('Distribución por mes')
ax.set_ylabel('Ventas (millones CLP)')

plt.tight_layout()
plt.show()
# Patrón estacional: ventas promedio por mes
df_mensual = pd.DataFrame({'ventas': serie.values, 'mes': serie.index.month,
                            'año': serie.index.year})
meses_nombre = ['Ene','Feb','Mar','Abr','May','Jun',
                'Jul','Ago','Sep','Oct','Nov','Dic']

fig, axes = plt.subplots(1, 2, figsize=(14, 4))

# Promedio por mes
ax = axes[0]
prom_mes = df_mensual.groupby('mes')['ventas'].mean()
colores_mes = ['crimson' if m == 12 else 'steelblue' for m in range(1, 13)]
ax.bar(range(1, 13), prom_mes/1e6, color=colores_mes, edgecolor='white', alpha=0.85)
ax.set_xticks(range(1, 13))
ax.set_xticklabels(meses_nombre)
ax.set_title('Ventas promedio por mes (2017–2023)')
ax.set_ylabel('Ventas (millones CLP)')

# Boxplot por mes
ax = axes[1]
data_box = [df_mensual[df_mensual['mes'] == m]['ventas'].values / 1e6
            for m in range(1, 13)]
bp = ax.boxplot(data_box, patch_artist=True, notch=False)
for patch, color in zip(bp['boxes'], colores_mes):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)
ax.set_xticklabels(meses_nombre)
ax.set_title('Distribución por mes')
ax.set_ylabel('Ventas (millones CLP)')

plt.tight_layout()
plt.show()

3. Descomposición de la Serie¶

Una serie temporal puede descomponerse en tres componentes:

$$Y_t = T_t + S_t + R_t \quad \text{(modelo aditivo)}$$ $$Y_t = T_t \times S_t \times R_t \quad \text{(modelo multiplicativo)}$$

Componente	Descripción
Tendencia ($T_t$)	Movimiento de largo plazo
Estacionalidad ($S_t$)	Patrón cíclico de período fijo
Residuo ($R_t$)	Lo que queda después de extraer tendencia y estacionalidad

💡 Usar modelo multiplicativo cuando la amplitud de la estacionalidad crece con el nivel de la serie. Aditivo cuando la amplitud es constante.

In [30]:

Copied!





descomp = seasonal_decompose(serie, model='additive', period=12)

fig, axes = plt.subplots(4, 1, figsize=(13, 10), sharex=True)

componentes = [
    (serie,              'Serie original',  'steelblue'),
    (descomp.trend,      'Tendencia',       'seagreen'),
    (descomp.seasonal,   'Estacionalidad',  'darkorange'),
    (descomp.resid,      'Residuo',         'gray'),
]

for ax, (datos, titulo, color) in zip(axes, componentes):
    ax.plot(datos.index, datos/1e6, color=color, lw=1.8)
    ax.set_ylabel('M CLP')
    ax.set_title(titulo, fontsize=10)
    if titulo == 'Residuo':
        ax.axhline(0, color='black', lw=0.8, linestyle='--')

axes[-1].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
plt.suptitle('Descomposición aditiva de la serie de ventas', y=1.01, fontsize=12)
plt.tight_layout()
plt.show()
descomp = seasonal_decompose(serie, model='additive', period=12)

fig, axes = plt.subplots(4, 1, figsize=(13, 10), sharex=True)

componentes = [
    (serie,              'Serie original',  'steelblue'),
    (descomp.trend,      'Tendencia',       'seagreen'),
    (descomp.seasonal,   'Estacionalidad',  'darkorange'),
    (descomp.resid,      'Residuo',         'gray'),
]

for ax, (datos, titulo, color) in zip(axes, componentes):
    ax.plot(datos.index, datos/1e6, color=color, lw=1.8)
    ax.set_ylabel('M CLP')
    ax.set_title(titulo, fontsize=10)
    if titulo == 'Residuo':
        ax.axhline(0, color='black', lw=0.8, linestyle='--')

axes[-1].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
plt.suptitle('Descomposición aditiva de la serie de ventas', y=1.01, fontsize=12)
plt.tight_layout()
plt.show()

4. Estacionariedad y Prueba ADF¶

La mayoría de los modelos de series temporales requieren que la serie sea estacionaria: media, varianza y autocorrelación constantes en el tiempo. Una serie con tendencia NO es estacionaria.

La Prueba de Dickey-Fuller Aumentada (ADF) contrasta:

$H_0$: la serie tiene raíz unitaria (no es estacionaria)
$H_1$: la serie es estacionaria

Si p-value < 0.05 → rechazamos $H_0$ → la serie es estacionaria.

In [31]:

Copied!





def test_adf(serie, nombre='Serie'):
    resultado = adfuller(serie.dropna(), autolag='AIC')
    print(f"=== Prueba ADF: {nombre} ===")
    print(f"  Estadístico ADF: {resultado[0]:.4f}")
    print(f"  p-value:         {resultado[1]:.4f}")
    print(f"  Valores críticos: ", {k: f"{v:.3f}" for k, v in resultado[4].items()})
    if resultado[1] < 0.05:
        print("  → ✅ Estacionaria (p < 0.05)")
    else:
        print("  → ⚠️  No estacionaria (p ≥ 0.05) — necesita diferenciación")
    print()

test_adf(serie, 'Ventas originales')

# Primera diferencia: elimina la tendencia
serie_diff = serie.diff().dropna()
test_adf(serie_diff, 'Primera diferencia')
def test_adf(serie, nombre='Serie'):
    resultado = adfuller(serie.dropna(), autolag='AIC')
    print(f"=== Prueba ADF: {nombre} ===")
    print(f"  Estadístico ADF: {resultado[0]:.4f}")
    print(f"  p-value:         {resultado[1]:.4f}")
    print(f"  Valores críticos: ", {k: f"{v:.3f}" for k, v in resultado[4].items()})
    if resultado[1] < 0.05:
        print("  → ✅ Estacionaria (p < 0.05)")
    else:
        print("  → ⚠️  No estacionaria (p ≥ 0.05) — necesita diferenciación")
    print()

test_adf(serie, 'Ventas originales')

# Primera diferencia: elimina la tendencia
serie_diff = serie.diff().dropna()
test_adf(serie_diff, 'Primera diferencia')

=== Prueba ADF: Ventas originales ===
  Estadístico ADF: 0.6174
  p-value:         0.9880
  Valores críticos:  {'1%': '-3.525', '5%': '-2.903', '10%': '-2.589'}
  → ⚠️  No estacionaria (p ≥ 0.05) — necesita diferenciación

=== Prueba ADF: Primera diferencia ===
  Estadístico ADF: -6.2924
  p-value:         0.0000
  Valores críticos:  {'1%': '-3.527', '5%': '-2.904', '10%': '-2.589'}
  → ✅ Estacionaria (p < 0.05)

In [32]:

Copied!





# Visualización: serie original vs diferenciada
fig, axes = plt.subplots(2, 1, figsize=(13, 6))

axes[0].plot(serie.index, serie/1e6, color='steelblue', lw=1.8)
axes[0].set_title('Serie original (no estacionaria)')
axes[0].set_ylabel('Ventas (M CLP)')

axes[1].plot(serie_diff.index, serie_diff/1e6, color='seagreen', lw=1.5)
axes[1].axhline(0, color='black', lw=0.8, linestyle='--')
axes[1].set_title('Primera diferencia (estacionaria)')
axes[1].set_ylabel('Δ Ventas (M CLP)')
axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

plt.tight_layout()
plt.show()
# Visualización: serie original vs diferenciada
fig, axes = plt.subplots(2, 1, figsize=(13, 6))

axes[0].plot(serie.index, serie/1e6, color='steelblue', lw=1.8)
axes[0].set_title('Serie original (no estacionaria)')
axes[0].set_ylabel('Ventas (M CLP)')

axes[1].plot(serie_diff.index, serie_diff/1e6, color='seagreen', lw=1.5)
axes[1].axhline(0, color='black', lw=0.8, linestyle='--')
axes[1].set_title('Primera diferencia (estacionaria)')
axes[1].set_ylabel('Δ Ventas (M CLP)')
axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

plt.tight_layout()
plt.show()

5. ACF y PACF¶

Los gráficos de Autocorrelación (ACF) y Autocorrelación Parcial (PACF) muestran cómo se relaciona la serie con sus valores pasados. Son la herramienta principal para identificar los órdenes del modelo ARIMA.

Gráfico	¿Qué muestra?	Cómo leerlo para ARIMA
ACF	Correlación entre $Y_t$ e $Y_{t-k}$ para cada lag $k$	Orden MA ($q$): se corta en el lag $q$
PACF	Correlación directa entre $Y_t$ e $Y_{t-k}$ (eliminando intermedios)	Orden AR ($p$): se corta en el lag $p$

In [33]:

Copied!





fig, axes = plt.subplots(2, 2, figsize=(14, 7))

# ACF y PACF de la serie original
plot_acf(serie,       lags=36, ax=axes[0,0], title='ACF — Serie original',       alpha=0.05)
plot_pacf(serie,      lags=36, ax=axes[0,1], title='PACF — Serie original',      alpha=0.05)

# ACF y PACF de la primera diferencia
plot_acf(serie_diff,  lags=36, ax=axes[1,0], title='ACF — Primera diferencia',   alpha=0.05)
plot_pacf(serie_diff, lags=36, ax=axes[1,1], title='PACF — Primera diferencia',  alpha=0.05)

for ax in axes.ravel():
    ax.set_xlabel('Lag (meses)')

plt.suptitle('ACF y PACF — Diagnóstico de autocorrelación', y=1.01, fontsize=12)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(2, 2, figsize=(14, 7))

# ACF y PACF de la serie original
plot_acf(serie,       lags=36, ax=axes[0,0], title='ACF — Serie original',       alpha=0.05)
plot_pacf(serie,      lags=36, ax=axes[0,1], title='PACF — Serie original',      alpha=0.05)

# ACF y PACF de la primera diferencia
plot_acf(serie_diff,  lags=36, ax=axes[1,0], title='ACF — Primera diferencia',   alpha=0.05)
plot_pacf(serie_diff, lags=36, ax=axes[1,1], title='PACF — Primera diferencia',  alpha=0.05)

for ax in axes.ravel():
    ax.set_xlabel('Lag (meses)')

plt.suptitle('ACF y PACF — Diagnóstico de autocorrelación', y=1.01, fontsize=12)
plt.tight_layout()
plt.show()

Qué observar:

En la serie original: picos significativos en lags 12, 24, 36 → estacionalidad anual clara
En la primera diferencia: la estructura se simplifica → la diferenciación eliminó la tendencia
Los picos en lag 12 de la diferencia → necesitamos diferenciación estacional también

6. Modelo SARIMA¶

El modelo SARIMA es la extensión de ARIMA para series con estacionalidad:

$$\text{SARIMA}(p,d,q)(P,D,Q)_s$$

Parámetro	Significado
$p$	Orden autorregresivo (lags AR)
$d$	Grado de diferenciación
$q$	Orden de media móvil
$P, D, Q$	Componentes estacionales equivalentes
$s$	Período estacional (12 para datos mensuales)

💡 En la práctica no es necesario calcular $(p,d,q)$ a mano. Lo importante es entender qué hace cada componente y saber leer el diagnóstico del modelo.

In [34]:

Copied!





# Split temporal: los últimos 12 meses como test
n_test = 12
serie_train = serie[:-n_test]
serie_test  = serie[-n_test:]

print(f"Train: {serie_train.index[0].date()} a {serie_train.index[-1].date()} ({len(serie_train)} obs)")
print(f"Test:  {serie_test.index[0].date()} a {serie_test.index[-1].date()}  ({len(serie_test)} obs)")
# Split temporal: los últimos 12 meses como test
n_test = 12
serie_train = serie[:-n_test]
serie_test  = serie[-n_test:]

print(f"Train: {serie_train.index[0].date()} a {serie_train.index[-1].date()} ({len(serie_train)} obs)")
print(f"Test:  {serie_test.index[0].date()} a {serie_test.index[-1].date()}  ({len(serie_test)} obs)")

Train: 2017-01-01 a 2022-12-01 (72 obs)
Test:  2023-01-01 a 2023-12-01  (12 obs)

In [35]:

Copied!





# Ajuste SARIMA(1,1,1)(1,1,1)[12]
modelo_sarima = SARIMAX(
    serie_train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

print(modelo_sarima.summary().tables[1])
# Ajuste SARIMA(1,1,1)(1,1,1)[12]
modelo_sarima = SARIMAX(
    serie_train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

print(modelo_sarima.summary().tables[1])

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1287      0.418     -0.308      0.758      -0.948       0.691
ma.L1         -0.7227      0.179     -4.030      0.000      -1.074      -0.371
ar.S.L12      -0.2907      0.291     -0.999      0.318      -0.861       0.279
ma.S.L12      -0.1763      0.274     -0.643      0.520      -0.714       0.361
sigma2      1.299e+10   1.02e-11   1.27e+21      0.000     1.3e+10     1.3e+10
==============================================================================

In [36]:

Copied!





# Diagnóstico de residuos del modelo
modelo_sarima.plot_diagnostics(figsize=(13, 7))
plt.suptitle('Diagnóstico de residuos — SARIMA(1,1,1)(1,1,1)[12]', y=1.01)
plt.tight_layout()
plt.show()
# Diagnóstico de residuos del modelo
modelo_sarima.plot_diagnostics(figsize=(13, 7))
plt.suptitle('Diagnóstico de residuos — SARIMA(1,1,1)(1,1,1)[12]', y=1.01)
plt.tight_layout()
plt.show()

In [37]:

Copied!





# Pronóstico sobre el período test
forecast_sarima = modelo_sarima.get_forecast(steps=n_test)
pred_mean = forecast_sarima.predicted_mean
pred_ci   = forecast_sarima.conf_int(alpha=0.05)

# Métricas
mae_s  = mean_absolute_error(serie_test, pred_mean)
rmse_s = mean_squared_error(serie_test, pred_mean) ** 0.5  # ← corregido
mape_s = (np.abs((serie_test - pred_mean) / serie_test)).mean() * 100

print("=== SARIMA — Métricas en test ===")
print(f"MAE:  {mae_s:,.0f} CLP")
print(f"RMSE: {rmse_s:,.0f} CLP")
print(f"MAPE: {mape_s:.2f}%")
# Pronóstico sobre el período test
forecast_sarima = modelo_sarima.get_forecast(steps=n_test)
pred_mean = forecast_sarima.predicted_mean
pred_ci   = forecast_sarima.conf_int(alpha=0.05)

# Métricas
mae_s  = mean_absolute_error(serie_test, pred_mean)
rmse_s = mean_squared_error(serie_test, pred_mean) ** 0.5  # ← corregido
mape_s = (np.abs((serie_test - pred_mean) / serie_test)).mean() * 100

print("=== SARIMA — Métricas en test ===")
print(f"MAE:  {mae_s:,.0f} CLP")
print(f"RMSE: {rmse_s:,.0f} CLP")
print(f"MAPE: {mape_s:.2f}%")

=== SARIMA — Métricas en test ===
MAE:  107,279 CLP
RMSE: 131,154 CLP
MAPE: 7.71%

In [38]:

Copied!





# Visualización del pronóstico SARIMA
fig, ax = plt.subplots(figsize=(13, 5))

ax.plot(serie_train.index, serie_train/1e6,
        color='steelblue', lw=1.8, label='Train')
ax.plot(serie_test.index, serie_test/1e6,
        color='black', lw=2, label='Real (test)')
ax.plot(pred_mean.index, pred_mean/1e6,
        color='crimson', lw=2, linestyle='--', label=f'SARIMA (MAPE={mape_s:.1f}%)')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0]/1e6,
                pred_ci.iloc[:, 1]/1e6,
                alpha=0.2, color='crimson', label='IC 95%')
ax.axvline(serie_test.index[0], color='gray', lw=1.5, linestyle='--')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.set_title('Pronóstico SARIMA — Últimos 12 meses')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
# Visualización del pronóstico SARIMA
fig, ax = plt.subplots(figsize=(13, 5))

ax.plot(serie_train.index, serie_train/1e6,
        color='steelblue', lw=1.8, label='Train')
ax.plot(serie_test.index, serie_test/1e6,
        color='black', lw=2, label='Real (test)')
ax.plot(pred_mean.index, pred_mean/1e6,
        color='crimson', lw=2, linestyle='--', label=f'SARIMA (MAPE={mape_s:.1f}%)')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0]/1e6,
                pred_ci.iloc[:, 1]/1e6,
                alpha=0.2, color='crimson', label='IC 95%')
ax.axvline(serie_test.index[0], color='gray', lw=1.5, linestyle='--')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.set_title('Pronóstico SARIMA — Últimos 12 meses')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

7. Enfoque ML: Features de Calendario y Lags¶

Una alternativa moderna a SARIMA es tratar el problema de forecasting como un problema de regresión supervisada. La clave está en construir features que capturen la información temporal:

Tipo de feature	Ejemplos	Captura
Features de calendario	mes, trimestre, día_semana, es_diciembre	Estacionalidad
Lags	$Y_{t-1}, Y_{t-2}, \ldots, Y_{t-12}$	Dependencia temporal
Estadísticas rodantes	media_3m, std_6m, max_12m	Tendencia y volatilidad reciente
Features de tendencia	$t$, $t^2$	Crecimiento a largo plazo

💡 La ventaja del enfoque ML es que puede incorporar fácilmente variables externas (precio del petróleo, indicadores económicos, campañas de marketing) que un SARIMA no maneja bien.

In [39]:

Copied!





def construir_features(serie, n_lags=12, ventanas=[3, 6, 12]):
    """
    Construye un DataFrame con features temporales para forecasting con ML.
    """
    df = pd.DataFrame({'y': serie})

    # Features de calendario
    df['mes']        = df.index.month
    df['trimestre']  = df.index.quarter
    df['es_dic']     = (df.index.month == 12).astype(int)
    df['es_enero']   = (df.index.month == 1).astype(int)
    df['tendencia']  = np.arange(len(df))

    # Lags
    for lag in range(1, n_lags + 1):
        df[f'lag_{lag}'] = df['y'].shift(lag)

    # Estadísticas rodantes
    for v in ventanas:
        df[f'media_{v}m']  = df['y'].shift(1).rolling(v).mean()
        df[f'std_{v}m']    = df['y'].shift(1).rolling(v).std()

    return df.dropna()


df_features = construir_features(serie, n_lags=12, ventanas=[3, 6, 12])

print(f"Features construidas: {df_features.shape[1] - 1}")
print(f"Observaciones disponibles: {len(df_features)}")
print(df_features.columns.tolist())
def construir_features(serie, n_lags=12, ventanas=[3, 6, 12]):
    """
    Construye un DataFrame con features temporales para forecasting con ML.
    """
    df = pd.DataFrame({'y': serie})

    # Features de calendario
    df['mes']        = df.index.month
    df['trimestre']  = df.index.quarter
    df['es_dic']     = (df.index.month == 12).astype(int)
    df['es_enero']   = (df.index.month == 1).astype(int)
    df['tendencia']  = np.arange(len(df))

    # Lags
    for lag in range(1, n_lags + 1):
        df[f'lag_{lag}'] = df['y'].shift(lag)

    # Estadísticas rodantes
    for v in ventanas:
        df[f'media_{v}m']  = df['y'].shift(1).rolling(v).mean()
        df[f'std_{v}m']    = df['y'].shift(1).rolling(v).std()

    return df.dropna()


df_features = construir_features(serie, n_lags=12, ventanas=[3, 6, 12])

print(f"Features construidas: {df_features.shape[1] - 1}")
print(f"Observaciones disponibles: {len(df_features)}")
print(df_features.columns.tolist())

Features construidas: 23
Observaciones disponibles: 72
['y', 'mes', 'trimestre', 'es_dic', 'es_enero', 'tendencia', 'lag_1', 'lag_2', 'lag_3', 'lag_4', 'lag_5', 'lag_6', 'lag_7', 'lag_8', 'lag_9', 'lag_10', 'lag_11', 'lag_12', 'media_3m', 'std_3m', 'media_6m', 'std_6m', 'media_12m', 'std_12m']

In [40]:

Copied!





# Split temporal: los últimos 12 meses como test
# IMPORTANTE: en series temporales el split debe respetar el orden cronológico
feature_cols = [c for c in df_features.columns if c != 'y']
X_all = df_features[feature_cols]
y_all = df_features['y']

X_train_ml = X_all.iloc[:-n_test]
X_test_ml  = X_all.iloc[-n_test:]
y_train_ml = y_all.iloc[:-n_test]
y_test_ml  = y_all.iloc[-n_test:]

print(f"Train ML: {len(X_train_ml)} observaciones")
print(f"Test ML:  {len(X_test_ml)} observaciones")
# Split temporal: los últimos 12 meses como test
# IMPORTANTE: en series temporales el split debe respetar el orden cronológico
feature_cols = [c for c in df_features.columns if c != 'y']
X_all = df_features[feature_cols]
y_all = df_features['y']

X_train_ml = X_all.iloc[:-n_test]
X_test_ml  = X_all.iloc[-n_test:]
y_train_ml = y_all.iloc[:-n_test]
y_test_ml  = y_all.iloc[-n_test:]

print(f"Train ML: {len(X_train_ml)} observaciones")
print(f"Test ML:  {len(X_test_ml)} observaciones")

Train ML: 60 observaciones
Test ML:  12 observaciones

In [41]:

Copied!





# Dos modelos ML para comparar
modelos_ml = {
    'Ridge + features':         Pipeline([('s', StandardScaler()), ('m', Ridge(alpha=100))]),
    'Gradient Boosting':        GradientBoostingRegressor(
                                    n_estimators=200, max_depth=3,
                                    learning_rate=0.05, random_state=42),
}

resultados_ml = {}
for nombre, modelo in modelos_ml.items():
    modelo.fit(X_train_ml, y_train_ml)
    pred = modelo.predict(X_test_ml)
    mae  = mean_absolute_error(y_test_ml, pred)
    rmse = mean_squared_error(y_test_ml, pred) ** 0.5  # ← corregido
    mape = (np.abs((y_test_ml - pred) / y_test_ml)).mean() * 100
    resultados_ml[nombre] = {"pred": pred, "MAE": mae, "RMSE": rmse, "MAPE": mape}
    print(f"[{nombre}]  MAE={mae:,.0f}  RMSE={rmse:,.0f}  MAPE={mape:.2f}%")
# Dos modelos ML para comparar
modelos_ml = {
    'Ridge + features':         Pipeline([('s', StandardScaler()), ('m', Ridge(alpha=100))]),
    'Gradient Boosting':        GradientBoostingRegressor(
                                    n_estimators=200, max_depth=3,
                                    learning_rate=0.05, random_state=42),
}

resultados_ml = {}
for nombre, modelo in modelos_ml.items():
    modelo.fit(X_train_ml, y_train_ml)
    pred = modelo.predict(X_test_ml)
    mae  = mean_absolute_error(y_test_ml, pred)
    rmse = mean_squared_error(y_test_ml, pred) ** 0.5  # ← corregido
    mape = (np.abs((y_test_ml - pred) / y_test_ml)).mean() * 100
    resultados_ml[nombre] = {"pred": pred, "MAE": mae, "RMSE": rmse, "MAPE": mape}
    print(f"[{nombre}]  MAE={mae:,.0f}  RMSE={rmse:,.0f}  MAPE={mape:.2f}%")

[Ridge + features]  MAE=71,684  RMSE=89,820  MAPE=5.32%
[Gradient Boosting]  MAE=78,106  RMSE=99,563  MAPE=5.72%

In [42]:

Copied!





# Importancia de features — Gradient Boosting
gb_model = modelos_ml['Gradient Boosting']
importancias = pd.Series(
    gb_model.feature_importances_,
    index=feature_cols
).sort_values(ascending=True).tail(15)

fig, ax = plt.subplots(figsize=(9, 5))
importancias.plot(kind='barh', color='steelblue', edgecolor='white', ax=ax)
ax.set_title('Top 15 features más importantes — Gradient Boosting')
ax.set_xlabel('Importancia relativa')
plt.tight_layout()
plt.show()
# Importancia de features — Gradient Boosting
gb_model = modelos_ml['Gradient Boosting']
importancias = pd.Series(
    gb_model.feature_importances_,
    index=feature_cols
).sort_values(ascending=True).tail(15)

fig, ax = plt.subplots(figsize=(9, 5))
importancias.plot(kind='barh', color='steelblue', edgecolor='white', ax=ax)
ax.set_title('Top 15 features más importantes — Gradient Boosting')
ax.set_xlabel('Importancia relativa')
plt.tight_layout()
plt.show()

8. Comparación Final: SARIMA vs ML¶

Evaluamos todos los modelos sobre el mismo período de test.

In [43]:

Copied!





# Tabla comparativa
filas = [{'Modelo': 'SARIMA(1,1,1)(1,1,1)[12]',
          'MAE': mae_s, 'RMSE': rmse_s, 'MAPE (%)': mape_s}]
for nombre, res in resultados_ml.items():
    filas.append({'Modelo': nombre,
                  'MAE': res['MAE'], 'RMSE': res['RMSE'], 'MAPE (%)': res['MAPE']})

df_comp = pd.DataFrame(filas).set_index('Modelo')
df_comp['MAE']  = df_comp['MAE'].map('{:,.0f}'.format)
df_comp['RMSE'] = df_comp['RMSE'].map('{:,.0f}'.format)
df_comp['MAPE (%)'] = df_comp['MAPE (%)'].round(2)
print("=== Comparación de modelos — Test set (últimos 12 meses) ===")
print(df_comp.to_string())
# Tabla comparativa
filas = [{'Modelo': 'SARIMA(1,1,1)(1,1,1)[12]',
          'MAE': mae_s, 'RMSE': rmse_s, 'MAPE (%)': mape_s}]
for nombre, res in resultados_ml.items():
    filas.append({'Modelo': nombre,
                  'MAE': res['MAE'], 'RMSE': res['RMSE'], 'MAPE (%)': res['MAPE']})

df_comp = pd.DataFrame(filas).set_index('Modelo')
df_comp['MAE']  = df_comp['MAE'].map('{:,.0f}'.format)
df_comp['RMSE'] = df_comp['RMSE'].map('{:,.0f}'.format)
df_comp['MAPE (%)'] = df_comp['MAPE (%)'].round(2)
print("=== Comparación de modelos — Test set (últimos 12 meses) ===")
print(df_comp.to_string())

=== Comparación de modelos — Test set (últimos 12 meses) ===
                              MAE     RMSE  MAPE (%)
Modelo                                              
SARIMA(1,1,1)(1,1,1)[12]  107,279  131,154      7.71
Ridge + features           71,684   89,820      5.32
Gradient Boosting          78,106   99,563      5.72

In [44]:

Copied!





# Visualización comparativa de pronósticos
fig, ax = plt.subplots(figsize=(13, 5))

# Últimos 24 meses de train para contexto
ax.plot(serie_train.index[-24:], serie_train.iloc[-24:]/1e6,
        color='steelblue', lw=1.8, label='Histórico')
ax.plot(serie_test.index, serie_test/1e6,
        color='black', lw=2.5, label='Real')

# Pronósticos
ax.plot(pred_mean.index, pred_mean/1e6,
        color='crimson', lw=2, linestyle='--',
        label=f"SARIMA (MAPE={mape_s:.1f}%)")

estilos = ['-.', ':']
colores_ml = ['seagreen', 'darkorange']
for (nombre, res), ls, color in zip(resultados_ml.items(), estilos, colores_ml):
    ax.plot(serie_test.index, res['pred']/1e6,
            color=color, lw=2, linestyle=ls,
            label=f"{nombre} (MAPE={res['MAPE']:.1f}%)")

ax.axvline(serie_test.index[0], color='gray', lw=1.5, linestyle='--')
ax.text(serie_test.index[0], ax.get_ylim()[0]*1.02, ' → Test',
        fontsize=9, color='gray')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.set_title('Comparación de pronósticos — SARIMA vs ML')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
# Visualización comparativa de pronósticos
fig, ax = plt.subplots(figsize=(13, 5))

# Últimos 24 meses de train para contexto
ax.plot(serie_train.index[-24:], serie_train.iloc[-24:]/1e6,
        color='steelblue', lw=1.8, label='Histórico')
ax.plot(serie_test.index, serie_test/1e6,
        color='black', lw=2.5, label='Real')

# Pronósticos
ax.plot(pred_mean.index, pred_mean/1e6,
        color='crimson', lw=2, linestyle='--',
        label=f"SARIMA (MAPE={mape_s:.1f}%)")

estilos = ['-.', ':']
colores_ml = ['seagreen', 'darkorange']
for (nombre, res), ls, color in zip(resultados_ml.items(), estilos, colores_ml):
    ax.plot(serie_test.index, res['pred']/1e6,
            color=color, lw=2, linestyle=ls,
            label=f"{nombre} (MAPE={res['MAPE']:.1f}%)")

ax.axvline(serie_test.index[0], color='gray', lw=1.5, linestyle='--')
ax.text(serie_test.index[0], ax.get_ylim()[0]*1.02, ' → Test',
        fontsize=9, color='gray')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.set_title('Comparación de pronósticos — SARIMA vs ML')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

9. Pronóstico hacia el Futuro¶

Finalmente, generamos un pronóstico para los próximos 12 meses usando el mejor modelo, entrenado sobre toda la serie disponible.

In [45]:

Copied!





# Re-entrenamos SARIMA con TODA la serie
modelo_final = SARIMAX(
    serie,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

horizonte = 12
forecast_futuro = modelo_final.get_forecast(steps=horizonte)
pred_futuro = forecast_futuro.predicted_mean
ci_futuro   = forecast_futuro.conf_int(alpha=0.05)

# Visualización del pronóstico futuro
fig, ax = plt.subplots(figsize=(13, 5))

ax.plot(serie.index, serie/1e6, color='steelblue', lw=1.8, label='Histórico')
ax.plot(pred_futuro.index, pred_futuro/1e6,
        color='crimson', lw=2.5, linestyle='--', label='Pronóstico 2024')
ax.fill_between(ci_futuro.index,
                ci_futuro.iloc[:, 0]/1e6,
                ci_futuro.iloc[:, 1]/1e6,
                alpha=0.2, color='crimson', label='IC 95%')

ax.axvline(pred_futuro.index[0], color='gray', lw=1.5, linestyle='--')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax.set_title('Pronóstico 2024 — SARIMA entrenado en serie completa')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

print("=== Pronóstico mensual 2024 ===")
resumen_pron = pd.DataFrame({
    'Predicción': pred_futuro,
    'IC inf':     ci_futuro.iloc[:, 0],
    'IC sup':     ci_futuro.iloc[:, 1]
})
resumen_pron.index = resumen_pron.index.strftime('%Y-%m')
print((resumen_pron / 1e6).round(2).to_string())
# Re-entrenamos SARIMA con TODA la serie
modelo_final = SARIMAX(
    serie,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

horizonte = 12
forecast_futuro = modelo_final.get_forecast(steps=horizonte)
pred_futuro = forecast_futuro.predicted_mean
ci_futuro   = forecast_futuro.conf_int(alpha=0.05)

# Visualización del pronóstico futuro
fig, ax = plt.subplots(figsize=(13, 5))

ax.plot(serie.index, serie/1e6, color='steelblue', lw=1.8, label='Histórico')
ax.plot(pred_futuro.index, pred_futuro/1e6,
        color='crimson', lw=2.5, linestyle='--', label='Pronóstico 2024')
ax.fill_between(ci_futuro.index,
                ci_futuro.iloc[:, 0]/1e6,
                ci_futuro.iloc[:, 1]/1e6,
                alpha=0.2, color='crimson', label='IC 95%')

ax.axvline(pred_futuro.index[0], color='gray', lw=1.5, linestyle='--')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax.set_title('Pronóstico 2024 — SARIMA entrenado en serie completa')
ax.set_xlabel('Fecha')
ax.set_ylabel('Ventas (millones CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

print("=== Pronóstico mensual 2024 ===")
resumen_pron = pd.DataFrame({
    'Predicción': pred_futuro,
    'IC inf':     ci_futuro.iloc[:, 0],
    'IC sup':     ci_futuro.iloc[:, 1]
})
resumen_pron.index = resumen_pron.index.strftime('%Y-%m')
print((resumen_pron / 1e6).round(2).to_string())

=== Pronóstico mensual 2024 ===
         Predicción  IC inf  IC sup
2024-01        1.26    1.03    1.50
2024-02        1.34    1.10    1.57
2024-03        1.21    0.97    1.45
2024-04        1.48    1.24    1.72
2024-05        1.60    1.36    1.85
2024-06        1.69    1.43    1.94
2024-07        1.71    1.45    1.96
2024-08        1.62    1.35    1.88
2024-09        1.61    1.34    1.87
2024-10        1.52    1.25    1.79
2024-11        1.51    1.24    1.79
2024-12        1.69    1.41    1.96

Resumen¶

Concepto	Idea central
Serie temporal	Secuencia ordenada en el tiempo; el orden importa para la modelación
Descomposición	Separar tendencia, estacionalidad y residuo para entender la estructura
Estacionariedad	La mayoría de modelos la requieren; se verifica con la prueba ADF
Diferenciación	Transforma una serie no estacionaria en estacionaria
ACF / PACF	Herramientas para diagnosticar autocorrelación e identificar órdenes ARIMA
SARIMA	Modelo clásico para series con tendencia y estacionalidad; interpretable
Enfoque ML	Construir features de calendario, lags y estadísticas rodantes; flexible y extensible
MAPE	Métrica de error en porcentaje; fácil de comunicar a audiencias no técnicas
Split temporal	Siempre cronológico: train primero, test al final (nunca aleatorio)

¿SARIMA o ML?¶

Situación	Recomendación
Serie univariada, estacionalidad clara, horizonte corto	SARIMA
Variables externas disponibles (precio, clima, marketing)	ML con features
Muchas series que modelar simultáneamente	ML (más escalable)
Necesitas intervalos de confianza interpretables	SARIMA
Patrones no lineales o múltiples estacionalidades	ML (Gradient Boosting)

Referencias¶

Hyndman, R.J., Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3ª ed.). https://otexts.com/fpp3/ — referencia gratuita en línea.
Statsmodels — SARIMAX: https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html
Scikit-learn — GradientBoostingRegressor: https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html