Exemple complet: regressió lineal pas a pas

La millor manera d’entendre un mètode és aplicar-lo de cap a cua sobre dades reals. En aquest subcapítol treballarem amb una sèrie de vendes mensuals sintètica i construirem un model de regressió lineal des de l’exploració inicial fins al diagnòstic final de residus.

Cada pas inclou el codi complet i l’explicació del perquè, no només del què.

Les dades: 60 mesos de vendes

Farem servir una sèrie de vendes mensuals de 5 anys amb tres components: tendència creixent, estacionalitat anual i soroll aleatori.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(42)
dates = pd.date_range(start='2020-01-01', periods=60, freq='ME')
t = np.arange(60)

tendència = 50 + 0.8 * t
estacionalitat = 10 * np.sin(2 * np.pi * t / 12)
soroll = np.random.normal(0, 4, 60)

vendes = tendència + estacionalitat + soroll
serie = pd.Series(vendes, index=dates, name='vendes')

Pas 1: Exploració inicial

print(f'Períodes: {len(serie)}')
print(f'Rang temporal: {serie.index[0].date()} → {serie.index[-1].date()}')
print(f'Mitjana: {serie.mean():.1f}')
print(f'Desv. típica: {serie.std():.1f}')

creixement_anual = (serie.iloc[-1] - serie.iloc[0]) / (len(serie) / 12)
print(f'Creixement anual estimat: {creixement_anual:.1f} unitats')

Pas 2: Construir el dataset de features

df = pd.DataFrame({'vendes': serie})
df['t'] = np.arange(len(df))

df['lag_1'] = df['vendes'].shift(1)
df['lag_2'] = df['vendes'].shift(2)
df['lag_12'] = df['vendes'].shift(12)
df['mitj_mov_3'] = df['vendes'].shift(1).rolling(3).mean()

df['mes'] = df.index.month
dummies_mes = pd.get_dummies(df['mes'], prefix='mes', drop_first=True)
df = pd.concat([df, dummies_mes], axis=1)

df_model = df.dropna()
print(f'Files disponibles per al model: {len(df_model)}')

Pas 3: Split temporal train/test

HORITZÓ = 12
train_size = len(df_model) - HORITZÓ

feature_cols = ['t', 'lag_1', 'lag_2', 'lag_12', 'mitj_mov_3'] + \
               [c for c in df_model.columns if c.startswith('mes_')]

X = df_model[feature_cols]
y = df_model['vendes']

X_train = X.iloc[:train_size]
X_test  = X.iloc[train_size:]
y_train = y.iloc[:train_size]
y_test  = y.iloc[train_size:]

Pas 4: Ajustar el model

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_train, y_train)

y_pred_train = model.predict(X_train)
y_pred_test  = model.predict(X_test)

coefs = pd.Series(model.coef_, index=feature_cols)
print('\nCoeficients del model:')
print(coefs.sort_values(key=abs, ascending=False).round(3))

Pas 5: Avaluar les prediccions

from sklearn.metrics import mean_absolute_error, mean_squared_error

mae_train = mean_absolute_error(y_train, y_pred_train)
mae_test  = mean_absolute_error(y_test, y_pred_test)
rmse_test = mean_squared_error(y_test, y_pred_test) ** 0.5

baseline = pd.Series(y_train.iloc[-1], index=y_test.index)
mae_baseline = mean_absolute_error(y_test, baseline)

print(f'MAE train:    {mae_train:.2f}')
print(f'MAE test:     {mae_test:.2f}')
print(f'RMSE test:    {rmse_test:.2f}')
print(f'MAE baseline: {mae_baseline:.2f}')
print(f'\nMillora sobre baseline: {(1 - mae_test/mae_baseline)*100:.1f}%')

Pas 6: Diagnosticar els residus

from statsmodels.stats.stattools import durbin_watson

residus_train = y_train - y_pred_train
dw = durbin_watson(residus_train)
print(f'Durbin-Watson: {dw:.3f}')

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(residus_train.index, residus_train, color='#1a1a1a', linewidth=1.2)
axes[0].axhline(0, color='#ffcc00', linestyle='--', linewidth=2)
axes[0].set_title(f'Residus en el temps (DW = {dw:.2f})')
axes[0].grid(True, alpha=0.3)

axes[1].hist(residus_train, bins=15, color='#ffcc00', edgecolor='#1a1a1a')
axes[1].set_title('Distribució de residus')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Pas 7: Millorar el model

df['lag_3'] = df['vendes'].shift(3)
df['lag_6'] = df['vendes'].shift(6)
df['t2'] = df['t'] ** 2

df_model_v2 = df.dropna()
feature_cols_v2 = ['t', 't2', 'lag_1', 'lag_2', 'lag_3', 'lag_6', 'lag_12',
                   'mitj_mov_3'] + [c for c in df_model_v2.columns if c.startswith('mes_')]

X_v2 = df_model_v2[feature_cols_v2]
y_v2 = df_model_v2['vendes']

train_size_v2 = len(df_model_v2) - HORITZÓ
model_v2 = LinearRegression()
model_v2.fit(X_v2.iloc[:train_size_v2], y_v2.iloc[:train_size_v2])

y_pred_v2 = model_v2.predict(X_v2.iloc[train_size_v2:])
mae_v2 = mean_absolute_error(y_v2.iloc[train_size_v2:], y_pred_v2)

print(f'MAE model v1: {mae_test:.2f}')
print(f'MAE model v2: {mae_v2:.2f}')

Resultats esperats

Què hem après

En aquest exemple has recorregut el cicle complet d’un model de regressió per a sèries temporals:

1. Explorar la sèrie visualment abans de modelar
2. Construir features temporals: lags, mitjanes mòbils, dummies estacionals
3. Separar train/test respectant l’ordre temporal
4. Ajustar un model de regressió lineal amb sklearn
5. Avaluar amb MAE, RMSE i comparació contra baseline
6. Diagnosticar residus amb Durbin-Watson i visualització
7. Iterar afegint features i validant en test

Aquest flux de treball és pràcticament idèntic al que faràs servir amb models més sofisticats com Random Forest o XGBoost per a sèries temporals. La diferència és en el model, no en el pipeline.