Analisi delle Serie Temporali in Python: Guida Pratica da pandas al Forecasting Moderno

Introduzione

Se lavori con i dati, prima o poi ti imbatterai nelle serie temporali. È praticamente inevitabile. Una serie temporale non è altro che una sequenza di osservazioni ordinate nel tempo: i prezzi giornalieri di un'azione in borsa, le temperature orarie di una stazione meteo, le vendite mensili di un prodotto, il consumo energetico di una città ogni quindici minuti. Ovunque ci sia un fenomeno che evolve nel tempo e viene misurato a intervalli regolari (o anche irregolari), ecco una serie temporale.

Ma perché dovresti preoccupartene? Beh, per almeno due ottime ragioni.

In primo luogo, analizzare una serie temporale ti permette di capire cosa sta succedendo davvero sotto la superficie: c'è un trend crescente? Una stagionalità che si ripete ogni settimana o ogni anno? In secondo luogo — e qui viene il bello dal punto di vista pratico — puoi fare previsioni sul futuro. Il cosiddetto forecasting è cruciale in finanza, logistica, energia, marketing e in decine di altri ambiti.

Python si è affermato come il linguaggio di riferimento per questo tipo di analisi, grazie a un ecosistema di librerie davvero impressionante. In questa guida percorreremo l'intero percorso: si parte dai fondamenti con pandas, si passa alla visualizzazione con matplotlib, alla decomposizione e ai test di stazionarietà con statsmodels, per poi arrivare ai modelli ARIMA e SARIMAX. Infine, esploreremo le soluzioni moderne dell'ecosistema Nixtla — StatsForecast e NeuralForecast — e chiuderemo con un progetto pratico end-to-end di previsione della domanda energetica.

Allacciati la cintura, sarà un viaggio lungo ma ne vale la pena.

Fondamenti delle Serie Temporali con pandas

La libreria pandas è il punto di partenza naturale per qualsiasi lavoro con serie temporali in Python. Il suo DataFrame, combinato con il potente sistema di indicizzazione temporale, offre strumenti davvero sofisticati per caricare, manipolare e trasformare dati temporali.

DatetimeIndex e parsing delle date

Il cuore della gestione delle serie temporali in pandas è il DatetimeIndex. La funzione pd.to_datetime() è lo strumento principale per convertire stringhe, numeri o altri formati in oggetti datetime — e onestamente, funziona sorprendentemente bene nella maggior parte dei casi.

import pandas as pd
import numpy as np

# Creazione di un DatetimeIndex da stringhe
date_strings = ['2024-01-01', '2024-01-02', '2024-01-03', '2024-01-04']
date_index = pd.to_datetime(date_strings)
print(date_index)

# Creazione di un intervallo di date con pd.date_range
date_range = pd.date_range(start='2024-01-01', periods=365, freq='D')
print(f"Intervallo: da {date_range[0]} a {date_range[-1]}")

# Creazione di una Serie Temporale sintetica
np.random.seed(42)
ts = pd.Series(
    data=np.cumsum(np.random.randn(365)) + 100,
    index=date_range,
    name='valore'
)
print(ts.head(10))

Quando carichi dati da un CSV, puoi specificare direttamente la colonna da usare come indice temporale durante la lettura. Questo ti evita conversioni successive ed è decisamente più pulito.

# Lettura di un CSV con parsing automatico delle date
df = pd.read_csv(
    'dati_vendite.csv',
    parse_dates=['data'],
    index_col='data'
)

# Parsing personalizzato per formati non standard
df['timestamp'] = pd.to_datetime(
    df['data_raw'],
    format='%d/%m/%Y %H:%M',
    dayfirst=True
)

# Gestione dei fusi orari
df.index = df.index.tz_localize('Europe/Rome')

Il parametro freq e la regolarità temporale

Il parametro freq definisce la frequenza della serie temporale ed è fondamentale per molte operazioni. Le frequenze più comuni includono 'D' (giornaliero), 'H' (orario), 'M' (fine mese), 'MS' (inizio mese), 'W' (settimanale), 'Q' (trimestrale) e 'Y' (annuale). Se l'indice non ha una frequenza esplicita, puoi impostarla o inferirla.

# Impostare la frequenza dell'indice
ts = ts.asfreq('D')

# Inferire la frequenza
inferred_freq = pd.infer_freq(ts.index)
print(f"Frequenza inferita: {inferred_freq}")

# Creare range con frequenze particolari
orario = pd.date_range('2024-01-01', periods=48, freq='h')
ogni_15_min = pd.date_range('2024-01-01', periods=96, freq='15min')
giorni_lavorativi = pd.date_range('2024-01-01', periods=60, freq='B')

Resampling: cambiare la frequenza temporale

Il resampling è una delle operazioni più potenti di pandas per le serie temporali. Ti consente di modificare la frequenza dei dati, sia aggregandoli (downsampling) sia interpolandoli (upsampling). Personalmente, è una delle funzionalità che uso più spesso nei miei progetti.

# Dati orari sintetici di consumo energetico
hourly_index = pd.date_range('2024-01-01', periods=8760, freq='h')
energy = pd.Series(
    np.random.uniform(200, 800, 8760) +
    100 * np.sin(np.arange(8760) * 2 * np.pi / 24),
    index=hourly_index,
    name='consumo_kwh'
)

# Downsampling: da orario a giornaliero
daily_mean = energy.resample('D').mean()
daily_sum = energy.resample('D').sum()
daily_max = energy.resample('D').max()

# Downsampling: da orario a mensile
monthly = energy.resample('MS').agg(['mean', 'sum', 'min', 'max'])
print(monthly.head())

# Upsampling: da giornaliero a orario con interpolazione
upsampled = daily_mean.resample('h').interpolate(method='cubic')

Rolling windows, shift e differencing

Le finestre mobili (rolling windows) sono fondamentali per smussare il rumore e identificare trend. L'operazione di shift crea variabili ritardate (lag), mentre il differencing è essenziale per rendere stazionaria una serie. Sono tre strumenti che, una volta imparati, userai continuamente.

# Media mobile a 7 e 30 giorni
ts_rolling_7 = ts.rolling(window=7).mean()
ts_rolling_30 = ts.rolling(window=30).mean()

# Media mobile esponenziale (EWM)
ts_ewm = ts.ewm(span=20).mean()

# Rolling con statistiche multiple
rolling_stats = ts.rolling(window=30).agg(['mean', 'std', 'min', 'max'])
print(rolling_stats.tail())

# Shift: creare lag per analisi di autocorrelazione
lag_1 = ts.shift(1)    # valore del giorno precedente
lag_7 = ts.shift(7)    # valore di una settimana fa
lag_365 = ts.shift(365) # valore di un anno fa

# Calcolo del rendimento percentuale
rendimento = ts.pct_change()

# Differencing di primo ordine
diff_1 = ts.diff(1)

# Differencing stagionale (es. stagionalità settimanale)
diff_seasonal = ts.diff(7)

# Differencing di secondo ordine
diff_2 = ts.diff(1).diff(1)

Visualizzazione delle Serie Temporali

La visualizzazione è un passaggio che non puoi permetterti di saltare. Ti permette di identificare a colpo d'occhio trend, stagionalità, anomalie e pattern che i numeri da soli non riescono a comunicare con la stessa immediatezza. Matplotlib, insieme a pandas, ti dà tutto quello che serve.

Grafici di base con matplotlib

import matplotlib.pyplot as plt
import matplotlib.dates as mdates

fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

# Grafico della serie originale con medie mobili
axes[0].plot(ts.index, ts, alpha=0.5, label='Dati originali', linewidth=0.8)
axes[0].plot(ts_rolling_7.index, ts_rolling_7, label='Media mobile 7gg', linewidth=2)
axes[0].plot(ts_rolling_30.index, ts_rolling_30, label='Media mobile 30gg', linewidth=2)
axes[0].set_title('Serie Temporale con Medie Mobili')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Grafico della distribuzione mensile (boxplot)
monthly_data = ts.groupby(ts.index.month)
axes[1].boxplot(
    [group.values for _, group in monthly_data],
    labels=['Gen','Feb','Mar','Apr','Mag','Giu',
            'Lug','Ago','Set','Ott','Nov','Dic']
)
axes[1].set_title('Distribuzione per Mese')
axes[1].grid(True, alpha=0.3)

# Grafico delle differenze prime
axes[2].plot(diff_1.index, diff_1, color='green', alpha=0.7, linewidth=0.8)
axes[2].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[2].set_title('Differenze Prime (Stazionarietà)')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('analisi_serie_temporale.png', dpi=150, bbox_inches='tight')
plt.show()

Grafici ACF e PACF

I grafici di autocorrelazione (ACF) e autocorrelazione parziale (PACF) sono strumenti diagnostici fondamentali per scegliere i parametri dei modelli ARIMA. L'ACF mostra la correlazione tra la serie e le sue versioni ritardate; il PACF, invece, mostra solo la correlazione diretta, rimuovendo l'effetto dei ritardi intermedi. Sembrano complessi, ma con un po' di pratica diventano abbastanza intuitivi.

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

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

# ACF: utile per identificare il parametro q (MA)
plot_acf(ts.dropna(), lags=40, ax=axes[0], alpha=0.05)
axes[0].set_title('Funzione di Autocorrelazione (ACF)')

# PACF: utile per identificare il parametro p (AR)
plot_pacf(ts.dropna(), lags=40, ax=axes[1], alpha=0.05, method='ywm')
axes[1].set_title('Funzione di Autocorrelazione Parziale (PACF)')

plt.tight_layout()
plt.show()

Decomposizione e Stazionarietà

Due concetti che devi assolutamente padroneggiare nell'analisi delle serie temporali sono la decomposizione e la stazionarietà. La decomposizione ti permette di scomporre una serie nelle sue componenti, mentre la stazionarietà è un requisito necessario per la maggior parte dei modelli statistici classici.

Decomposizione stagionale

Una serie temporale può essere vista come la combinazione di tre componenti: il trend (la tendenza di lungo periodo), la stagionalità (pattern che si ripetono a intervalli regolari) e il residuo (la componente casuale, quello che resta dopo aver tolto trend e stagionalità). La decomposizione può essere additiva (Y = T + S + R) o moltiplicativa (Y = T × S × R), a seconda della natura della serie.

from statsmodels.tsa.seasonal import seasonal_decompose, STL

# Creiamo una serie con trend e stagionalità evidenti
np.random.seed(42)
t = np.arange(365)
trend = 0.05 * t
seasonal = 10 * np.sin(2 * np.pi * t / 365) + 5 * np.sin(2 * np.pi * t / 7)
noise = np.random.normal(0, 2, 365)
ts_complex = pd.Series(
    trend + seasonal + noise + 50,
    index=pd.date_range('2024-01-01', periods=365, freq='D'),
    name='valore'
)

# Decomposizione additiva classica
decomposition = seasonal_decompose(ts_complex, model='additive', period=30)

fig, axes = plt.subplots(4, 1, figsize=(14, 10))
decomposition.observed.plot(ax=axes[0], title='Serie Osservata')
decomposition.trend.plot(ax=axes[1], title='Trend')
decomposition.seasonal.plot(ax=axes[2], title='Stagionalità')
decomposition.resid.plot(ax=axes[3], title='Residui')
for ax in axes:
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Decomposizione STL (più robusta e flessibile)
stl = STL(ts_complex, period=30, robust=True)
result_stl = stl.fit()
result_stl.plot()
plt.show()

Test di stazionarietà: Augmented Dickey-Fuller

Una serie è stazionaria quando le sue proprietà statistiche — media, varianza, autocorrelazione — non cambiano nel tempo. La maggior parte dei modelli classici lo richiede. Il test ADF (Augmented Dickey-Fuller) è quello più utilizzato: l'ipotesi nulla dice che la serie ha una radice unitaria (cioè non è stazionaria). Se il p-value è inferiore a 0.05, puoi rifiutare l'ipotesi nulla e concludere che la serie è stazionaria.

In parole povere: p-value basso, serie stazionaria. Facile, no?

from statsmodels.tsa.stattools import adfuller, kpss

def test_stazionarieta(series, nome='Serie'):
    """Esegue il test ADF e KPSS su una serie temporale."""
    print(f"
{'='*50}")
    print(f"Test di stazionarietà per: {nome}")
    print(f"{'='*50}")

    # Test Augmented Dickey-Fuller
    result_adf = adfuller(series.dropna(), autolag='AIC')
    print(f"
Test ADF:")
    print(f"  Statistica del test: {result_adf[0]:.4f}")
    print(f"  P-value: {result_adf[1]:.6f}")
    print(f"  Lag utilizzati: {result_adf[2]}")
    print(f"  Valori critici:")
    for key, value in result_adf[4].items():
        print(f"    {key}: {value:.4f}")

    if result_adf[1] < 0.05:
        print(f"  => La serie È stazionaria (p < 0.05)")
    else:
        print(f"  => La serie NON è stazionaria (p >= 0.05)")

    return result_adf[1]

# Test sulla serie originale
p_orig = test_stazionarieta(ts_complex, 'Serie originale')

# Test dopo differencing di primo ordine
p_diff1 = test_stazionarieta(ts_complex.diff(1), 'Dopo diff(1)')

# Test dopo differencing di secondo ordine
p_diff2 = test_stazionarieta(ts_complex.diff(1).diff(1), 'Dopo diff(2)')

Il differencing è il metodo standard per rendere stazionaria una serie. Il numero di differenze necessarie corrisponde al parametro d nei modelli ARIMA. Nella pratica, una o due differenze sono quasi sempre sufficienti — raramente ne servono di più. Se la serie presenta stagionalità, potrebbe essere necessario anche un differencing stagionale.

Modelli Classici con statsmodels

Statsmodels è la libreria di riferimento in Python per i modelli statistici classici delle serie temporali. I modelli ARIMA e la loro estensione stagionale SARIMAX rappresentano il fondamento della modellazione statistica. Sono in giro da decenni, ma rimangono incredibilmente efficaci.

Il modello ARIMA

Un modello ARIMA è definito da tre parametri: p (ordine della componente autoregressiva), d (ordine di differenziazione) e q (ordine della media mobile). La selezione di questi parametri può essere guidata dall'analisi ACF/PACF e da criteri informativi come AIC e BIC.

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
warnings.filterwarnings('ignore')

# Preparazione dei dati: serie mensile di vendite
np.random.seed(42)
months = pd.date_range('2020-01-01', periods=48, freq='MS')
trend = np.linspace(100, 180, 48)
seasonal = 20 * np.sin(np.arange(48) * 2 * np.pi / 12)
noise = np.random.normal(0, 5, 48)
vendite = pd.Series(trend + seasonal + noise, index=months, name='vendite')

# Divisione train/test
train = vendite[:'2023-06']
test = vendite['2023-07':]
print(f"Training: {len(train)} mesi, Test: {len(test)} mesi")

# Fitting di un modello ARIMA(1,1,1)
model_arima = ARIMA(train, order=(1, 1, 1))
result_arima = model_arima.fit()
print(result_arima.summary())

# Previsioni
forecast_arima = result_arima.forecast(steps=len(test))
print("
Previsioni ARIMA:")
print(forecast_arima)

Il modello SARIMAX

Quando la serie ha una componente stagionale evidente, il modello SARIMAX è la scelta giusta. Oltre ai parametri (p,d,q), include quelli stagionali (P,D,Q,s), dove s è il periodo stagionale (12 per dati mensili con stagionalità annuale, 7 per dati giornalieri con pattern settimanale, e così via). SARIMAX supporta anche variabili esogene, cioè predittori esterni che possono migliorare le previsioni.

# Fitting di un modello SARIMAX con stagionalità annuale
model_sarimax = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
)
result_sarimax = model_sarimax.fit(disp=False)
print(result_sarimax.summary())

# Previsioni con intervalli di confidenza
forecast_obj = result_sarimax.get_forecast(steps=len(test))
forecast_mean = forecast_obj.predicted_mean
forecast_ci = forecast_obj.conf_int(alpha=0.05)

# Visualizzazione delle previsioni
fig, ax = plt.subplots(figsize=(14, 6))
train.plot(ax=ax, label='Training', linewidth=1.5)
test.plot(ax=ax, label='Test (valori reali)', linewidth=1.5)
forecast_mean.plot(ax=ax, label='Previsioni SARIMAX',
                   linewidth=2, linestyle='--')
ax.fill_between(
    forecast_ci.index,
    forecast_ci.iloc[:, 0],
    forecast_ci.iloc[:, 1],
    alpha=0.2, color='orange',
    label='Intervallo di confidenza 95%'
)
ax.set_title('Previsione delle Vendite con SARIMAX')
ax.set_ylabel('Vendite')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Diagnostica del modello

Dopo aver adattato un modello, devi verificare che i residui si comportino bene: dovrebbero essere approssimativamente normali, non autocorrelati e con media vicina a zero. Se non è così, il modello probabilmente non sta catturando qualche pattern importante nei dati.

# Diagnostica dei residui
result_sarimax.plot_diagnostics(figsize=(14, 10))
plt.tight_layout()
plt.show()

# Selezione automatica dei parametri tramite grid search
import itertools

p = d = q = range(0, 3)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in pdq]

best_aic = float('inf')
best_params = None
best_seasonal = None

for param in pdq:
    for seasonal_param in seasonal_pdq:
        try:
            mod = SARIMAX(
                train,
                order=param,
                seasonal_order=seasonal_param,
                enforce_stationarity=False,
                enforce_invertibility=False
            )
            results = mod.fit(disp=False, maxiter=100)
            if results.aic < best_aic:
                best_aic = results.aic
                best_params = param
                best_seasonal = seasonal_param
        except Exception:
            continue

print(f"Miglior modello: SARIMAX{best_params}x{best_seasonal}")
print(f"AIC: {best_aic:.2f}")

Metriche di valutazione

Per valutare quantitativamente le previsioni, si usano metriche standard: MAE (Mean Absolute Error), RMSE (Root Mean Squared Error) e MAPE (Mean Absolute Percentage Error). Il MAPE è particolarmente intuitivo perché esprime l'errore in percentuale.

from sklearn.metrics import mean_absolute_error, mean_squared_error

def calcola_metriche(reali, previsti, nome_modello='Modello'):
    """Calcola e stampa le metriche di valutazione."""
    mae = mean_absolute_error(reali, previsti)
    rmse = np.sqrt(mean_squared_error(reali, previsti))
    mape = np.mean(np.abs((reali - previsti) / reali)) * 100

    print(f"
Metriche per {nome_modello}:")
    print(f"  MAE:  {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.2f}%")

    return {'MAE': mae, 'RMSE': rmse, 'MAPE': mape}

metriche_sarimax = calcola_metriche(test, forecast_mean, 'SARIMAX')

L'Ecosistema Nixtla: StatsForecast per il Forecasting su Larga Scala

Fin qui tutto bene, ma cosa succede quando passi da una singola serie temporale a centinaia — o migliaia — di serie da prevedere simultaneamente? Le implementazioni tradizionali di statsmodels diventano lente. Molto lente.

Ed è qui che entra in gioco StatsForecast, la libreria sviluppata da Nixtla. È specificamente progettata per il forecasting su larga scala e offre implementazioni dei modelli classici che sono fino a 20 volte più veloci rispetto a pmdarima. Non sto esagerando.

Installazione e struttura dei dati

StatsForecast usa una convenzione specifica per i dati: un DataFrame con tre colonne — unique_id per identificare ciascuna serie, ds per il timestamp e y per il valore osservato. Questo formato standardizzato permette di gestire naturalmente più serie nello stesso DataFrame.

# pip install statsforecast
from statsforecast import StatsForecast
from statsforecast.models import (
    AutoARIMA, AutoETS, AutoTheta, AutoCES,
    SeasonalNaive, Naive, HistoricAverage
)

# Preparazione dati nel formato richiesto da StatsForecast
np.random.seed(42)

records = []
for store_id in ['negozio_A', 'negozio_B', 'negozio_C']:
    dates = pd.date_range('2021-01-01', periods=156, freq='W-MON')
    base = np.random.uniform(80, 120)
    trend = np.linspace(0, 30, 156)
    seasonal = 15 * np.sin(np.arange(156) * 2 * np.pi / 52)
    noise = np.random.normal(0, 5, 156)
    values = base + trend + seasonal + noise

    for d, v in zip(dates, values):
        records.append({
            'unique_id': store_id,
            'ds': d,
            'y': max(v, 0)
        })

df_nixtla = pd.DataFrame(records)
print(df_nixtla.head(10))
print(f"
Serie uniche: {df_nixtla['unique_id'].nunique()}")
print(f"Osservazioni totali: {len(df_nixtla)}")

Modelli disponibili e API unificata

StatsForecast mette a disposizione una vasta gamma di modelli, ciascuno con una versione "Auto" che seleziona automaticamente i parametri ottimali. Ecco i principali:

• AutoARIMA: seleziona automaticamente i parametri (p,d,q)(P,D,Q,s) — equivalente a pmdarima, ma molto più veloce
• AutoETS: modelli di livellamento esponenziale con selezione automatica
• AutoTheta: implementazione del metodo Theta, particolarmente efficace per serie con trend
• AutoCES: Complex Exponential Smoothing, una variante moderna del livellamento esponenziale
• SeasonalNaive e Naive: modelli baseline semplici ma essenziali per il confronto (non sottovalutarli mai!)

# Definizione dei modelli da confrontare
models = [
    AutoARIMA(season_length=52),
    AutoETS(season_length=52),
    AutoTheta(season_length=52),
    AutoCES(season_length=52),
    SeasonalNaive(season_length=52),
]

# Creazione dell'oggetto StatsForecast
sf = StatsForecast(
    models=models,
    freq='W-MON',
    n_jobs=-1,        # utilizza tutti i core disponibili
    fallback_model=SeasonalNaive(season_length=52),
)

# Fitting e previsione in un unico passaggio
# horizon = 12 settimane di previsione
forecasts = sf.forecast(
    df=df_nixtla,
    h=12,
    level=[90, 95],   # intervalli di confidenza
)

print(forecasts.head(15))
print(f"
Colonne generate: {forecasts.columns.tolist()}")

Cross-validation con StatsForecast

Una delle funzionalità che apprezzo di più in StatsForecast è la cross-validation temporale integrata. Permette di valutare le prestazioni dei modelli su più finestre temporali mantenendo il corretto ordinamento cronologico — tutto con poche righe di codice.

# Cross-validation con finestre scorrevoli
cv_results = sf.cross_validation(
    df=df_nixtla,
    h=12,             # orizzonte di previsione
    step_size=4,      # spostamento tra finestre
    n_windows=5,      # numero di finestre di validazione
    level=[90],
)

print(cv_results.head())

# Calcolo delle metriche per ogni modello
model_names = ['AutoARIMA', 'AutoETS', 'AutoTheta', 'AutoCES', 'SeasonalNaive']

for model_name in model_names:
    if model_name in cv_results.columns:
        errors = cv_results['y'] - cv_results[model_name]
        mae = np.mean(np.abs(errors))
        rmse = np.sqrt(np.mean(errors**2))
        print(f"{model_name:20s} - MAE: {mae:.2f}, RMSE: {rmse:.2f}")

NeuralForecast: Modelli Deep Learning per le Serie Temporali

I modelli classici sono spesso competitivi (a volte sorprendentemente), ma quando i pattern diventano complessi, le relazioni non lineari e i dati abbondanti, il deep learning può fare la differenza. NeuralForecast, sempre dell'ecosistema Nixtla, fornisce implementazioni dei principali modelli neurali mantenendo la stessa interfaccia e formato dati di StatsForecast. Cosa che, detto tra noi, rende il passaggio dall'uno all'altro praticamente indolore.

NBEATS e LSTM

N-BEATS (Neural Basis Expansion Analysis for Time Series) è un'architettura neurale progettata specificamente per il forecasting, con risultati eccellenti in numerose competizioni. LSTM (Long Short-Term Memory) è una rete neurale ricorrente adatta a catturare dipendenze temporali a lungo termine.

# pip install neuralforecast
from neuralforecast import NeuralForecast
from neuralforecast.models import NBEATS, LSTM, NHITS, PatchTST

# Definizione dei modelli neurali
horizon = 12

neural_models = [
    NBEATS(
        input_size=2 * horizon,   # finestra di input
        h=horizon,                 # orizzonte di previsione
        max_steps=300,            # epoche di addestramento
        scaler_type='standard',
        accelerator='cpu',
    ),
    NHITS(
        input_size=2 * horizon,
        h=horizon,
        max_steps=300,
        scaler_type='standard',
        accelerator='cpu',
    ),
]

# Creazione e addestramento del modello
nf = NeuralForecast(
    models=neural_models,
    freq='W-MON',
)

# Addestramento e previsione
nf.fit(df=df_nixtla)
neural_forecasts = nf.predict()
print(neural_forecasts.head())

Quando usare il deep learning

La scelta tra modelli classici e deep learning non è sempre scontata. Ecco alcune linee guida pratiche che mi sono state utili:

• Preferisci i modelli classici quando hai poche serie, serie brevi (meno di qualche centinaio di osservazioni), pattern lineari e stagionalità regolare, o quando l'interpretabilità è fondamentale
• Considera il deep learning quando hai molte serie correlate, serie molto lunghe con pattern complessi, relazioni non lineari, o quando disponi di GPU
• Approccio consigliato: parti sempre da un baseline (Naive, SeasonalNaive) e da un modello classico (AutoARIMA, AutoETS). Solo dopo, se necessario, sperimenta con i modelli neurali per verificare se il guadagno giustifica la complessità aggiuntiva

Progetto Pratico: Previsione della Domanda Energetica

Bene, è arrivato il momento di mettere insieme tutto quello che abbiamo visto. Costruiremo un progetto end-to-end di previsione della domanda energetica per tre zone cittadine. Useremo un dataset sintetico ma realistico che simula il consumo orario, con pattern giornalieri e stagionali.

Generazione e preparazione dei dati

import pandas as pd
import numpy as np
from datetime import datetime

np.random.seed(42)

def genera_domanda_energetica(zona, n_giorni=730):
    """Genera dati sintetici di domanda energetica per una zona."""
    ore = pd.date_range('2022-01-01', periods=n_giorni*24, freq='h')
    n = len(ore)
    t = np.arange(n)

    # Trend lento crescente
    trend = 500 + 0.01 * t

    # Stagionalità annuale (estate più consumi per AC)
    giorno_anno = np.array([d.timetuple().tm_yday for d in ore])
    stag_annuale = 100 * np.sin(2 * np.pi * giorno_anno / 365 - np.pi/2)

    # Stagionalità giornaliera (picchi mattina e sera)
    ora = np.array([d.hour for d in ore])
    stag_giornaliera = (
        80 * np.sin(2 * np.pi * (ora - 6) / 24) +
        40 * np.sin(2 * np.pi * (ora - 18) / 12)
    )

    # Effetto weekend (meno consumo)
    giorno_settimana = np.array([d.weekday() for d in ore])
    effetto_weekend = np.where(giorno_settimana >= 5, -60, 0)

    # Rumore
    rumore = np.random.normal(0, 30, n)

    # Combinazione
    moltiplicatore = {'zona_nord': 1.0, 'zona_centro': 1.3, 'zona_sud': 0.85}
    domanda = (trend + stag_annuale + stag_giornaliera +
               effetto_weekend + rumore) * moltiplicatore.get(zona, 1.0)
    domanda = np.maximum(domanda, 50)

    return pd.DataFrame({
        'unique_id': zona,
        'ds': ore,
        'y': domanda
    })

# Generazione dati per tre zone
df_energia = pd.concat([
    genera_domanda_energetica('zona_nord'),
    genera_domanda_energetica('zona_centro'),
    genera_domanda_energetica('zona_sud'),
], ignore_index=True)

print(f"Shape del dataset: {df_energia.shape}")
print(f"Periodo: da {df_energia['ds'].min()} a {df_energia['ds'].max()}")
print(f"
Statistiche per zona:")
print(df_energia.groupby('unique_id')['y'].describe().round(1))

Analisi esplorativa

# Aggregazione giornaliera per visualizzazione
df_daily = df_energia.copy()
df_daily['data'] = df_daily['ds'].dt.date
daily_agg = df_daily.groupby(['unique_id', 'data'])['y'].mean().reset_index()

fig, axes = plt.subplots(3, 1, figsize=(16, 10), sharex=True)

for idx, zona in enumerate(['zona_nord', 'zona_centro', 'zona_sud']):
    dati_zona = daily_agg[daily_agg['unique_id'] == zona]
    axes[idx].plot(pd.to_datetime(dati_zona['data']),
                   dati_zona['y'], linewidth=0.7)
    axes[idx].set_title(f'Domanda Energetica Media Giornaliera - {zona}')
    axes[idx].set_ylabel('MW')
    axes[idx].grid(True, alpha=0.3)

plt.xlabel('Data')
plt.tight_layout()
plt.show()

# Analisi del profilo giornaliero medio
profilo_orario = df_energia.copy()
profilo_orario['ora'] = profilo_orario['ds'].dt.hour

fig, ax = plt.subplots(figsize=(12, 5))
for zona in df_energia['unique_id'].unique():
    dati = profilo_orario[profilo_orario['unique_id'] == zona]
    media_oraria = dati.groupby('ora')['y'].mean()
    ax.plot(media_oraria.index, media_oraria.values,
            marker='o', label=zona, linewidth=2)

ax.set_xlabel('Ora del giorno')
ax.set_ylabel('Domanda media (MW)')
ax.set_title('Profilo di Carico Giornaliero Medio per Zona')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xticks(range(24))
plt.tight_layout()
plt.show()

Resampling e decomposizione

# Resampling a frequenza giornaliera per i modelli
records_daily = []
for zona in ['zona_nord', 'zona_centro', 'zona_sud']:
    zona_data = df_energia[df_energia['unique_id'] == zona].copy()
    zona_data = zona_data.set_index('ds')['y'].resample('D').mean().reset_index()
    zona_data['unique_id'] = zona
    records_daily.append(zona_data)

df_model = pd.concat(records_daily, ignore_index=True)
df_model = df_model[['unique_id', 'ds', 'y']]

print(f"Dataset per il modelling: {df_model.shape}")
print(df_model.head())

# Decomposizione per una zona
zona_nord = df_model[df_model['unique_id'] == 'zona_nord'].set_index('ds')['y']
decomp = seasonal_decompose(zona_nord, model='additive', period=7)
decomp.plot()
plt.suptitle('Decomposizione - Zona Nord (periodo=7 giorni)')
plt.tight_layout()
plt.show()

Forecasting con StatsForecast

from statsforecast import StatsForecast
from statsforecast.models import (
    AutoARIMA, AutoETS, AutoTheta, SeasonalNaive
)

# Definizione dei modelli
horizon = 30  # previsione a 30 giorni

models_energy = [
    AutoARIMA(season_length=7),
    AutoETS(season_length=7),
    AutoTheta(season_length=7),
    SeasonalNaive(season_length=7),
]

sf_energy = StatsForecast(
    models=models_energy,
    freq='D',
    n_jobs=-1,
)

# Cross-validation per valutazione robusta
cv_energy = sf_energy.cross_validation(
    df=df_model,
    h=horizon,
    step_size=30,
    n_windows=3,
    level=[90],
)

# Calcolo metriche per ogni modello e zona
print("
Risultati Cross-Validation:")
print("=" * 65)

model_cols = ['AutoARIMA', 'AutoETS', 'AutoTheta', 'SeasonalNaive']
results_summary = []

for zona in df_model['unique_id'].unique():
    cv_zona = cv_energy[cv_energy['unique_id'] == zona]
    print(f"
{zona}:")
    for model_name in model_cols:
        if model_name in cv_zona.columns:
            errori = cv_zona['y'] - cv_zona[model_name]
            mae = np.mean(np.abs(errori))
            rmse = np.sqrt(np.mean(errori**2))
            mape = np.mean(np.abs(errori / cv_zona['y'])) * 100
            print(f"  {model_name:20s} MAE={mae:7.2f}  "
                  f"RMSE={rmse:7.2f}  MAPE={mape:5.2f}%")
            results_summary.append({
                'zona': zona, 'modello': model_name,
                'MAE': mae, 'RMSE': rmse, 'MAPE': mape
            })

df_results = pd.DataFrame(results_summary)

Previsione finale e visualizzazione

# Previsione finale sui dati completi
final_forecasts = sf_energy.forecast(
    df=df_model,
    h=horizon,
    level=[90],
)

# Visualizzazione delle previsioni per ogni zona
fig, axes = plt.subplots(3, 1, figsize=(16, 12))

for idx, zona in enumerate(['zona_nord', 'zona_centro', 'zona_sud']):
    ax = axes[idx]

    # Dati storici (ultimi 90 giorni)
    storico = df_model[df_model['unique_id'] == zona].tail(90)
    ax.plot(storico['ds'], storico['y'],
            label='Storico', color='black', linewidth=1)

    # Previsioni dei diversi modelli
    prev_zona = final_forecasts[final_forecasts['unique_id'] == zona]
    colors = {'AutoARIMA': 'blue', 'AutoETS': 'red',
              'AutoTheta': 'green', 'SeasonalNaive': 'orange'}

    for model_name, color in colors.items():
        if model_name in prev_zona.columns:
            ax.plot(prev_zona['ds'], prev_zona[model_name],
                    label=model_name, color=color,
                    linewidth=1.5, linestyle='--')

    # Intervallo di confidenza per AutoARIMA
    lo_col = 'AutoARIMA-lo-90'
    hi_col = 'AutoARIMA-hi-90'
    if lo_col in prev_zona.columns and hi_col in prev_zona.columns:
        ax.fill_between(
            prev_zona['ds'],
            prev_zona[lo_col],
            prev_zona[hi_col],
            alpha=0.15, color='blue'
        )

    ax.set_title(f'Previsione Domanda Energetica - {zona}')
    ax.set_ylabel('MW medio giornaliero')
    ax.legend(loc='upper left', fontsize=8)
    ax.grid(True, alpha=0.3)

plt.xlabel('Data')
plt.tight_layout()
plt.savefig('previsione_energia.png', dpi=150, bbox_inches='tight')
plt.show()

# Identificazione del miglior modello complessivo
best_per_zona = df_results.loc[df_results.groupby('zona')['MAE'].idxmin()]
print("
Miglior modello per zona (in base al MAE):")
for _, row in best_per_zona.iterrows():
    print(f"  {row['zona']}: {row['modello']} "
          f"(MAE={row['MAE']:.2f}, MAPE={row['MAPE']:.2f}%)")

Best Practice e Consigli

L'analisi delle serie temporali ha le sue insidie specifiche. Vediamo le best practice essenziali per non cadere nelle trappole più comuni (e credimi, ci sono cascato anch'io più di una volta).

Cross-validation per le serie temporali

A differenza della cross-validation classica, nelle serie temporali devi rispettare l'ordine cronologico. Non puoi mai usare dati futuri per prevedere il passato — sembra ovvio, ma è un errore più comune di quanto pensi. Scikit-learn fornisce TimeSeriesSplit, mentre StatsForecast ha la propria implementazione integrata.

from sklearn.model_selection import TimeSeriesSplit

# TimeSeriesSplit di scikit-learn per modelli custom
tscv = TimeSeriesSplit(n_splits=5, test_size=30)

serie_singola = df_model[
    df_model['unique_id'] == 'zona_nord'
].set_index('ds')['y']

fig, axes = plt.subplots(5, 1, figsize=(14, 10))

for fold, (train_idx, test_idx) in enumerate(tscv.split(serie_singola)):
    train_fold = serie_singola.iloc[train_idx]
    test_fold = serie_singola.iloc[test_idx]

    axes[fold].plot(train_fold.index, train_fold.values,
                    label='Train', color='blue')
    axes[fold].plot(test_fold.index, test_fold.values,
                    label='Test', color='red')
    axes[fold].set_title(f'Fold {fold + 1}: Train={len(train_fold)}, '
                        f'Test={len(test_fold)}')
    axes[fold].legend(loc='upper left', fontsize=8)
    axes[fold].grid(True, alpha=0.3)

plt.suptitle('TimeSeriesSplit - 5 Fold', fontsize=14)
plt.tight_layout()
plt.show()

Evitare il data leakage

Il data leakage è probabilmente il rischio più insidioso nell'analisi delle serie temporali. Si verifica quando informazioni dal futuro finiscono (involontariamente) nei dati di training. Ecco le regole d'oro:

• Scaling e normalizzazione vanno calcolati solo sui dati di training, poi applicati al test
• Le feature basate su statistiche aggregate devono usare solo dati fino al momento della previsione
• Mai usare split casuali: separa sempre cronologicamente training e test
• L'imputation dei valori mancanti deve rispettare l'ordine temporale (forward fill, non interpolazione bidirezionale)

# Esempio CORRETTO di preprocessing senza data leakage
from sklearn.preprocessing import StandardScaler

# Divisione cronologica
split_date = '2023-07-01'
train_data = serie_singola[serie_singola.index < split_date]
test_data = serie_singola[serie_singola.index >= split_date]

# Scaling calcolato SOLO sul training
scaler = StandardScaler()
train_scaled = scaler.fit_transform(train_data.values.reshape(-1, 1))
test_scaled = scaler.transform(test_data.values.reshape(-1, 1))

# Feature engineering con rolling SOLO su dati passati
def create_features_safe(data, window=7):
    """Crea feature basate su finestre mobili senza data leakage."""
    features = pd.DataFrame(index=data.index)
    features['lag_1'] = data.shift(1)
    features['lag_7'] = data.shift(7)
    features['rolling_mean_7'] = data.shift(1).rolling(window).mean()
    features['rolling_std_7'] = data.shift(1).rolling(window).std()
    # Nota: shift(1) assicura che usiamo solo dati passati
    return features

features_train = create_features_safe(train_data)
features_test = create_features_safe(
    pd.concat([train_data, test_data])  # servono i lag dal training
).loc[test_data.index]

Gestione dei valori mancanti

I valori mancanti sono la norma, non l'eccezione, nelle serie temporali reali. Guasti dei sensori, interruzioni di rete, errori di registrazione... le cause sono infinite. La scelta del metodo di imputazione dipende dalla natura dei dati e dalla lunghezza delle interruzioni.

# Simulazione di valori mancanti
serie_con_nan = serie_singola.copy()
np.random.seed(42)
idx_mancanti = np.random.choice(len(serie_con_nan), size=50, replace=False)
serie_con_nan.iloc[idx_mancanti] = np.nan
print(f"Valori mancanti: {serie_con_nan.isna().sum()}")

# Metodo 1: Forward fill (propagazione in avanti)
serie_ffill = serie_con_nan.ffill()

# Metodo 2: Interpolazione lineare
serie_interp = serie_con_nan.interpolate(method='linear')

# Metodo 3: Interpolazione temporale (tiene conto dell'indice)
serie_time_interp = serie_con_nan.interpolate(method='time')

# Metodo 4: Interpolazione con limite massimo di gap
serie_limited = serie_con_nan.interpolate(
    method='linear', limit=3, limit_direction='forward'
)

# Metodo 5: Sostituzione con valore stagionale
# (media dello stesso giorno della settimana)
giorno_sett = serie_con_nan.index.dayofweek
medie_giorno = serie_con_nan.groupby(giorno_sett).transform('mean')
serie_stagionale = serie_con_nan.fillna(medie_giorno)

# Verifica
print(f"
Valori mancanti dopo ffill: {serie_ffill.isna().sum()}")
print(f"Valori mancanti dopo interpolazione: {serie_interp.isna().sum()}")
print(f"Valori mancanti dopo media stagionale: {serie_stagionale.isna().sum()}")

Scegliere tra modelli

Come scegliere il modello giusto? Non c'è una risposta universale, ma ecco una guida pratica che funziona nella maggior parte dei casi:

1. Parti sempre da un baseline: Naive e SeasonalNaive sono il tuo punto di riferimento. Se un modello complesso non li supera in modo significativo, probabilmente non ne vale la pena
2. Considera la lunghezza della serie: per serie brevi (meno di 100 osservazioni), ETS o Theta sono generalmente preferibili. Per serie molto lunghe, modelli più complessi possono catturare pattern sottili
3. Valuta la stagionalità: se è forte e regolare, AutoETS e SeasonalNaive spesso bastano. Per stagionalità multiple o irregolari, AutoARIMA o modelli neurali possono essere più adatti
4. Numero di serie: per poche serie, statsmodels va benissimo. Per centinaia o migliaia, StatsForecast è la scelta naturale
5. Risorse computazionali: i modelli deep learning richiedono GPU e tempi di addestramento molto superiori. Assicurati che il guadagno in precisione giustifichi il costo
6. Interpretabilità: in contesti dove devi spiegare le previsioni (finanza regolamentata, medicina), i modelli classici sono preferibili per la loro trasparenza

Conclusione

In questa guida abbiamo attraversato l'intero spettro dell'analisi delle serie temporali in Python, dai fondamenti con pandas fino alle soluzioni moderne per il forecasting su larga scala. Facciamo un rapido riepilogo.

pandas resta lo strumento indispensabile per la manipolazione dei dati temporali. Indicizzazione temporale, resampling, rolling windows e differencing sono la base di tutto. Padroneggiare questi strumenti non è opzionale.

La visualizzazione non è un optional. Grafici della serie, decomposizioni e plot ACF/PACF ti danno informazioni che nessuna tabella di numeri può sostituire.

La comprensione della stazionarietà e della decomposizione è il prerequisito per applicare correttamente i modelli statistici. Il test ADF e il differencing devono far parte del tuo toolkit standard.

I modelli classici ARIMA e SARIMAX di statsmodels rimangono competitivi e affidabili. La loro solidità teorica e interpretabilità li rendono la prima scelta in molti contesti — non lasciarti sedurre troppo in fretta dal deep learning.

L'ecosistema Nixtla con StatsForecast è un salto di qualità per il forecasting operativo. Velocità, semplicità dell'API e supporto nativo per serie multiple lo rendono perfetto per la produzione.

NeuralForecast apre le porte al deep learning con architetture specializzate come N-BEATS e NHITS, capaci di catturare pattern che i modelli classici non possono cogliere — a patto di avere dati sufficienti e risorse computazionali adeguate.

Per chi vuole proseguire, vale la pena esplorare i modelli gerarchici (HierarchicalForecast di Nixtla), il rilevamento anomalie, l'integrazione in pipeline MLOps e i modelli foundation come TimeGPT, che rappresentano la nuova frontiera del settore.

L'analisi delle serie temporali è un campo in continua evoluzione, dove la teoria statistica classica si fonde con le innovazioni del machine learning. La chiave del successo? Scegliere lo strumento giusto per ciascun problema, partire sempre dalla comprensione dei dati, e validare rigorosamente ogni approccio. Buon forecasting!