Time Series: Analisi dei dati di traffico passeggeri di un aeroporto Italiano

12 min readDec 28, 2018

Rome Fiumicino Airport from the air (www.formiche.net)

Introduzione.

Le tecniche di analisi per le Serie Temporali possono essere applicate all’analisi ed alla predizione in molti campi dell’attività umana. Una delle aree tradizionalmente esplorata concerne i dati del traffico passeggeri degli aeroporti, ove è ragionevole ipotizzare la presenza di trend a medio-lungo termine e di variazioni dovute a stagionalità.

Tuttavia, in letteratura, si trovano soltanto analisi relative ad aeroporti stranieri (per me) e con dati non molto recenti. Per tale ragione ho provato a procurare e ad analizzare dati aggiornati relativi ad un aeroporto italiano: l’aeroporto di Fiumicino, il più importante aeroporto d’Italia e della città di Roma.

Alcune delle domande a cui volevo dare risposta con queste analisi:

qual è il trend del traffico passeggeri su Roma? E’ cambiato nel tempo? C’è crisi?
quali dati e quali strumenti di facile accesso sono disponibili per queste analisi?
quanto è difficile fare queste analisi?

Partiamo, come sempre, dai dati!

Il reperimento dei dati.

Il sito web di “Aeroporti di Roma (ADR)”, la società che gestisce gli scali aeroportuali romani, fornisce dati con cadenza mensile sia per l’aeroporto di Fiumicino (FCO) sia per l’aeroporto di Ciampino, che ha un traffico molto inferiore. I dati disponibili, accurati e dettagliati, sono relativi al periodo gennaio 2004 — novembre 2018.

Purtroppo, non è scaricabile un dataset direttamente utilizzabile e, pertanto, il primo compito che ho affrontato è come automatizzare l’acquisizione di tali dati, evitando la lunga e noiosa operazione manuale di accesso a più di un centinaio di pagine.

Per tale scopo, ho sviluppato un programma Python, che impiega una delle più usate librerie per fare “web scraping “: la libreria BeautifulSoap:

Con tale programma, il processo di invocare le URL corrispondenti ai 167 mesi (dal 2004 ad oggi), interpretare le pagine ed estrarre il dato mensile relativo al totale di passeggeri in transito, è completamente automatizzato e, re-direzionando l’output, genera un file in formato CSV, contenente i dati della Serie Temporale:

month, passenger“2004–01”,1785471“2004–02”,1835597“2004–03”,2203695“2004–04”,2354331“2004–05”,2461606“2004–06”,2481026“2004–07”,2769381“2004–08”,2790218“2004–09”,2699830....

Dopo aver prodotto il file, è indispensabile una prima analisi esplorativa allo scopo di controllare la qualità dei dati ed identificare eventuali anomalie. Ed infatti, il dato relativo ad aprile 2006 è risultato subito mancante. Ho potuto correggere l’anomalia leggendo manualmente la pagina web relativa e calcolando il numero di passeggeri dai dati parziali.

Il file CSV cosi’ prodotto e ripulito (“dati_fco2004.csv”) è la base di tutte le successive analisi.

Cosa possiamo procedere per analizzare la serie temporale?

L’obiettivo è di sviluppare un modello che consenta di interpretare quanto è avvenuto nel passato e possa essere utilizzato come punto di partenza per predizioni (forecast) future. In questo articolo mi occupo soltanto della parte di “Time Series Analysis”. Dedicherò, eventualmente, un secondo articolo allo sviluppo di modelli predittivi.

Quindi, seguendo un approccio classico, vogliamo sviluppare un modello che decomponga la serie temporale nelle seguenti componenti:

Trend
Stagionalità
Rumore (Noise)

Nello sviluppare il modello, si deve inoltre stabilire se è più opportuno adottare un modello additivo o moltiplicativo.

Modello additivo:

M(t) = T(t) + S(t) + N(t)

ove

M(t) è il generico valore, al tempo t, osservato
T(t) è la componente di trend
S(t) la componente di Stagionalità
N(t) il rumore

Modello moltiplicativo:

M(t) = T(t)*S(t)*N(t)

Strumenti utilizzati.

Il linguaggio di sviluppo utilizzato è Python, con l’impiego delle librerie Pandas, StatSmodels e Mathplotlib.

L’ambiente di sviluppo utilizzato è Jupyter, basato su Notebook. I vantaggi derivanti dall’utilizzo di un Notebook sono tanti:

E’ possibile combinare nello stesso file: codice, risultati delle elaborazioni (ad esempio, grafici e visualizzazioni di dati) e commenti per spiegare codice e risultati (markdown)
E’ agevole la ripetizione delle elaborazioni, dopo le modifiche suggerite dai risultati delle esplorazioni

L’approccio che ho adottato è quello tipico della moderna “Data Science”, un approccio esplorativo, basato su ipotesi, verifiche e rapide iterazioni per via via sviluppare e raffinare i modelli.

Seconda esplorazione.

Dopo la prima esplorazione del file, per eliminare dati che avrebbero creato problemi in caricamento, i dati del file CSV si possono caricare con un’unica istruzione in una Pandas Series, lo strumento di base per tutte le successive analisi.

FILE_NAME = “dati_fco2004.csv”series = read_csv(FILE_NAME, header=0, index_col = 0, parse_dates=True, squeeze=True)

La prima cosa che conviene sempre fare è produrre un grafico dei dati grezzi. Di seguito il risultato (sull’asse y il numero di passeggeri per mese).

Fig .1: Traffico passeggeri per mese (FCO 2004–2018)

Un grafico di questo tipo consente di:

escludere anomalie evidenti, importanti
capire se è presente un trend (è probabile, da quanto si vede)
evidenziare se è presente una variazione soggetta a stagionalità, periodica (anche questo è abbastanza evidente)

Un tale grafico, a mio, avviso, consente inoltre di decidere se è opportuno sviluppare un singolo modello, oppure se è opportuno dividere in sotto-serie e sviluppare modelli distinti. Dal grafico è abbastanza evidente che dopo un periodo di crescita (trend positivo) che dura fino ad inizio 2009, intervengono altre dinamiche che portano ad una crisi e successiva stasi nel periodo 2010–2013, per poi ritornare alla crescita. E’ opportuno, per avere maggiore accuratezza, segmentare i dati e sviluppare tre modelli distinti per i tre periodi.

Per confermare la presenza di una componente di stagionalità (seasonal) è utile produrre un grafico di auto-correlazione (auto-correlation plot):

Nell’”auto-correlation plot” si prendono in considerazione tutti i dati e si esamina la correlazione tra dati che abbiano lo stesso lag nel tempo (es: coppie di dati a 12 mesi di distanza, lag=12), al variare del lag. L’andamento sinusoidale osservabile nel grafico di sopra è chiaro indice di stagionalità e non sorprende che vi sia un elevata (coefficiente=0.8) correlazione positiva per lag=12. Le attività umane, e cosi i viaggi in aereo, seguono cicli a cadenza annuale, il massimo si ripete in estate e cosi via.

Analisi del primo periodo (2004–2009)

Nel primo periodo il trend di crescita è evidente ed univoco. L’ampiezza delle “oscillazioni” aumenta nel tempo e quindi non è efficace un modello additivo; Conviene usare un modello moltiplicativo:

M(t) = T(t)*S(t)*N(T)

Una decomposizione semplice e veloce si ottiene con la funzione seasonal_decompose() della libreria StatsModels

from statsmodels.tsa.seasonal import seasonal_decomposeresult2004 = seasonal_decompose(series2004, model=’multiplicative’)result2004.plot()

Fig . 4: Risultati della decomposizione (T,S,Residual)

I grafici confermano:

un trend di crescita abbastanza lineare, che però rallenta ad inizio 2008
una componente periodica, con periodo 12 (mesi)

Il rallentamento è spiegato dalla storia economica mondiale: è facile trovare su Internet i commenti relativi alla crisi economica ed alla crisi globale del traffico aereo iniziata nel 2008 e durata per tutto il 2009.

Un altro grafico interessante è quello che mostra il totale dei passeggeri in transito per anno. In esso si vede netto il calo del 2009 (ed anche quello del 2013).

Il valore, più basso, del 2018 si spiega con il fatto che nella serie temporale mancano i dati di dicembre 2018, non ancora disponibili.

Con un valore prevedibile intorno ai 3 milioni (passeggeri), ottenuto da un modello di forecasting SARIMA, si supera il dato del 2016 di circa 1 milione di passeggeri. Complimenti FCO!

Da quanto detto, è chiaro che un modello sviluppato sui dati 2004–2008 non avrebbe previsto correttamente il numero (in calo) di passeggeri del 2009. Ma, i modelli prevedono correttamente soltanto se non intervengono fenomeni ed eventi nuovi, non “contenuti” nei dati utilizzati per “l’addestramento” del modello stesso.

Il terzo periodo (2014–2018)

Proviamo anche in questo caso con un modello moltiplicativo.

Il risultato della decomposizione mediante seasonal_decompose() è rappresentato nei tre grafici seguenti:

Fig. 7: Risultati della decomposizione per il terzo periodo (T,S,Residual)

In questo caso l’analisi è differente. La componente di stagionalità ha sempre ampiezza e periodo costanti (è imposto dal metodo di decomposizione), ma il trend non è sempre crescente. Nel 2017 il trend si inverte, anche se l’effetto sul totale dei passeggeri è in parte compensato dalla componente di stagionalità (quest’ultima ha un massimo quando il trend raggiunge il minimo locale).

Cosa accade in quel periodo? Questa volta la crisi è “locale”. Il grafico mostra gli effetti della crisi della compagnia Alitalia, sopratutto nel traffico domestico, che colpisce principalmente il “suo hub” di Fiumicino.

Nel 2018 la crescita, trainata dagli effetti globali, riprende, almeno fino ad oggi.

Dal punto di vista della modellazione, per rappresentare un trend del tipo in figura, si deve ricorrere ad un “polynomial fitting”, con curva del terzo ordine.

Il secondo periodo (2010–2014)

Il grafico temporale mostra un periodo in cui non è evidente un trend di crescita e le variazioni sono in pratica legate alla sola componente stagionale.

Poiché l’ampiezza delle “oscillazioni” è praticamente costante, possiamo provare con un modello più semplice, un modello additivo:

M(t) = T(t) + S(t) + N(t)

il risultato della decomposizione è:

Fig. 9: Decomposizione (T,S,N), secondo periodo.

Questo risultato va attentamente analizzato, in quanto abbastanza problematico. E’, infatti, importante considerare:

il modello additivo impiegato
gli ordini di grandezza delle tre componenti (T, S, N) , in particolare il fatto che gli ordini di grandezza delle variazioni del trend e del rumore sono eguali (e quindi possono compensarsi)

In generale, sarei portato ad affermare che un modello in cui una componente di segnale (T) è dello stesso ordine di grandezza del rumore (N) non è molto utile. Nel caso specifico però l’equivalenza riguarda il trend e non anche la componente stagionale, molto più ampia. Il risultato è che in sostanza T + N prevedono un livello costante a meno delle fluttuazioni del rumore, a cui si sommano oscillazioni stagionali di ampiezza costante più ampia. E ciò coincide con quanto mostrato dal grafico della serie temporale del secondo periodo.

Approfondimento sulla componente stagionale del secondo periodo.

Laddove è presente una componente stagionale che produce alti e bassi, è spesso importante escluderla per analizzare con accuratezza il trend, per poter capire se nel medio e lungo periodo vi è una crescita complessiva. Una tecnica che si adotta è di modellare esplicitamente la componente stagionale e poi sottrarla alla serie originale.

Utilizzando i dati di questo periodo (2010–2014) applichiamo un’altra tecnica, utile per estrarre la componente stagionale e poter poi creare una serie temporale de-stagionalizzata (“seasonally adjusted”).

Avendo stabilito che la componente stagionale ha periodo 12, possiamo pensare di modellarla esplicitamente. Un comportamento periodico di ampiezza e periodo noti si può modellare con una funzione sinusoidale, che può essere adeguatamente approssimata con una funzione polinomiale di grado opportuno.

La libreria Numpy offre una funzione, polyfit() che può essere adattata ai dati (fit). Il grado che scegliamo per il polinomio è 4 (vedi “precisazioni”, dopo).

Il codice utilizzato per modellare la componente stagionale è il seguente:

Ed il grafico del modello polinomiale della componente stagionale è il seguente:

Fig .10: Secondo periodo (2010–2016): componente stagionale.

Per ottenere una serie temporale da cui è eliminata la componente stagionale basta fare la differenza:

series2010adj = series2010 - seasonal2010

Le frequenze della componente stagionale.

Proviamo ad applicare un ulteriore tecnica, tratta dalla Teoria dei Segnali.

L’idea è di trattare la serie temporale come un segnale somma di segnali sinusoidali a differente frequenza e di estrarre lo spettro del segnale effettuando una Trasformata Rapida di Fourier (FFT).

Per chiarezza:

prendiamo come unità temporale l’anno
la frequenza di campionamento usata è 12 (12 campioni/anno)
Una frequenza pari ad 1 vuol dire un segnale con un’oscillazione completa l’anno

Utilizzando la funzione fft del package scipy.fftpack si ottiene:

Fig. 11: FFT, secondo periodo (2010–2014)

Come interpretiamo il grafico?

Abbiamo una componente di frequenza 0 ed ampiezza di poco superiore a 3000000; Tale componente, costante, è il livello della serie
Poi, abbiamo una componente importante di frequenza 1: ciò conferma una significativa stagionalità con periodo pari ad 1 anno
Le altre componenti di frequenza 1/3 e 2/3 determinano una lenta variazione del trend, con periodi rispettivamente di 3 ed 1.5 anni.

Conclusioni.

Bene, una prima confortante conclusione (in questo periodo natalizio) è che il traffico passeggeri, dopo un periodo di stasi, ha ripreso a crescere, superando anche l’ultima crisi della compagnia di bandiera.

Ma, ovviamente, lo scopo di quest’articolo non era di scrivere una pagina inedita nella storia del traffico aereo italiano. Da questo punto di vista, con un pizzico di sana auto-ironia, potrei dire che, più che un articolo, si tratta di “un compitino”.

Tornando ad essere seri, l’obiettivo era di mostrare come tecniche di Time Series Analysis possono essere applicate, senza grosse difficoltà, a casi reali, avendo a disposizione i dati. Ed era di applicare tali analisi ad un dataset relativo ad un’importante realtà italiana.

Dal punto di vista delle tecniche adottate, in sintesi, abbiamo mostrato come:

estrarre i dati da un sito web, con un procedimento automatico
caricare i dati in una “Pandas Series”, utilizzando Python
esplorare i dati, producendo velocemente rappresentazioni grafiche
verificare la presenza di una componente periodica (seasonal), utilizzando un grafico di auto-correlazione
effettuare un’analisi spettrale della componente periodica
decomporre i dati nelle componenti di trend e stagionalità, per poterle analizzare separatamente
creare un modello esplicito (polinomiale) della componente stagionale

Alcune precisazioni.

Per i più precisi ed attenti, aggiungo alcuni particolari e giustificazioni.

I dati su cui ho lavorato sono salvati in un repository GitHub: https://github.com/luigisaetta/TimeSeriesfiumicino2018/blob/master/dati_fco2004.csv

Per quanto riguarda la modellazione polinomiale della componente stagionale, ho utilizzato un polinomio di grado 4. Perché 4 e non 3 o 5?

Il grado è un “iper-parametro” del modello. Per procedere in modo più rigoroso avrei dovuto applicare una “valutazione comparata su griglia”, utilizzando come metrica di accuratezza il RMSE. Per semplicità ho adottato un approccio più empirico.

In sintesi: in un periodo (12 mesi) dobbiamo avere un minimo ed un massimo. Poiché min, max sono gli zeri della derivata prima, ovviamente il polinomio deve avere come minimo grado 3. D’altro canto un polinomio di grado 4 ha 5 parametri. Più sono i parametri e più alto è il rischio di overfitting, tenuto conto del fatto che i dati non sono tanti (48). Il grado 4 adottato deriva da queste considerazioni e da un paio di tentativi con esplorazione visuale, applicando anche il principio del “keep it simple”.

Altra precisazione: molto spesso nell’effettuare la decomposizione nelle tre componenti (T, S, N), nei modelli moltiplicativi la componente S è normalizzata ad ampiezza 1 (vedi ad esempio il grafico di fig. 7). Ciò vuol dire che un fattore moltiplicativo è inglobato nel trend T, per dare il corretto ordine di grandezza complessivo. Nell’analisi della sezione dedicata ad approfondire la componente stagionale del secondo periodo si è scelto di mantenere in S il fattore e di non normalizzare.

Per quanto riguarda la FFT, si può notare dal grafico che non sono presenti frequenze superiori a 6. E’ conseguenza del teorema di Nyquist-Shannon, che stabilisce che: campionando ad un frequenza N non si possono ottenere informazioni su componenti di frequenza superiore ad N/2. (N = 12).

Per finire, per coloro i quali ci tengono ad essere sempre “sull’edge” della tecnologia: le tecniche utilizzate sono tecniche “classiche”, note in letteratura da parecchi anni. Oggi l’esplosione del “Deep Learning” porta ad analizzare le Serie Temporali con strumenti quali le reti LSTM. A riguardo credo che un opinione condivisibile sia di partire comunque da tecniche relativamente semplici, come quelle qui utilizzate, e poi di scegliere, se utile, algoritmi più complessi e computazionalmente più onerosi.

Previsione per dicembre 2018 (scritta il 29/12/2018).

I dati della serie storica, alla data in cui li ho raccolti (26 dicembre 2018) non contenevano ancora il dato relativo al mese di dicembre 2018.

Come potete immaginare, ho avuto la forte curiosità di provare ad applicare metodi di forecasting per prevedere tale valore, per poi confrontarlo con il dato reale quando disponibile.

Applicando un modello SARIMA, che tiene conto della stagionalità, la predizione che ho ottenuto è:

3011042 passeggeri (con un errore stimato in 78550)

Alcuni riferimenti per approfondire.