Una serie storica consta dei valori osservati in un insieme di periodi ordinati sequenzialmente. Questo, per chi fa SEO, è già un elemento del massimo interesse.
I dati di traffico del nostro sito web, considerati lungo una sequenza temporale, sono infatti un esempio di serie storica.
L’analisi delle serie storiche è un insieme di metodi che ci consentono di ricavare schemi o statistiche significative dai dati con informazioni temporali.
In termini molto generali, possiamo dire che una serie temporale è una sequenza di variabili casuali indicizzate nel tempo.
Lo scopo dell’analisi di una serie storica può essere di tipo descrittivo (si pensi alla decomposizione della serie per rimuovere elementi di stagionalità o per evidenziare tendenze di fondo) oppure inferenziale, includendo in quest’ultimo la previsione dei valori per periodi di tempo futuri, ancora non occorsi (forecasting).
Il metodo classico di analisi delle serie storiche individua quattro influenze, o componenti:
Secondo il modello dell’analisi classica delle serie storiche, il valore della variabile in ogni periodo è determinato dalle influenze delle quattro componenti.
Lo scopo principale dell’analisi classica delle serie temporali è proprio quello di scomporre la serie, per isolare le influenze delle varie componenti che determinano i valori della serie storica.
Le quattro componenti posso essere tra loro legate in modo additivo:
Y = T + C + S + I
ovvero in modo moltiplicativo:
Y = T x C x S x I
Ricordo che un modello moltiplicativo può essere trasformato nel modello additivo sfruttando le proprietà dei logaritmi:
log(Y) = log(T) + log(C) + log(S) + log(I)
Il logaritmo di un numero n nella base c (con c diverso da 1 e c >0) è l’esponente al quale è necessario elevare la base c per ottenere n.
Dunque, se n = cb allora logc n = b
Esistono vari modi per trasformare un vettore di dati, una matrice oppure un data frame in una serie temporale.
In questa sede ci limiteremo agli strumenti offerti, per ottenere tale risultato, dal pacchetto base R.
La funzione di nostro interesse si chiama semplicemente ts() ed ha un utilizzo piuttosto intuitivo.
Vediamo un esempio pratico. Supponiamo di avere un vettore di dati salvato con il nome mieidati.csv nella cartella /home del mio pc.
La prima cosa che farò, sarà quella di importare in R i miei dati che ipotizzo presenti in un file csv:
# importo i dati da file csv in un dataframe
dfmieidati <- read_csv(“/home/mieidati.csv”)
# creo un vettore con i dati che mi interessano
mieidati <- dfimieidati$miaosservazione
Ora non mi resta che invocare la funzione ts() con gli opportuni valori di inizio e di frequenza per creare la mia serie temporale.
Ipotizziamo che i dati siano mensili e comincino con gennaio 2012 e vadano fino a dicembre 2018:
serietemporale <- ts(mieidati, start=c(2012,1), end=c(2018,12), frequency=12)
come si vede la forma tipica della funzione ts() è
ts(vettore, inizio, fine, frequenza)
dove per frequenza si intende il numero di osservazioni per unità di tempo.
Quindi avremo 1= annuale, 4=trimestrale, 12=mensile…
Posso facilmente conoscere il tempo della prima osservazione di una serie storica usando il comando:
start()
Analogamente, posso trovare l’ultima osservazione con:
end()
Il comando:
frequency()
mi restituisce il numero delle osservazioni per unità di tempo, mentre l’utilissimo comando
window()
consente di estrarre un sottoinsieme di dati.
Usando a mo’ d’esempio un dataset compreso in R chiamato Nile contenente 100 letture annuali del fiume Nilo ad Aswan per gli anni dal 1871 al 197 (lo useremo anche nel paragrafo che segue) il comando si presenta così:
nile_sub <- window(Nile, start=1940,end=1960)
Uno dei vantaggi nell’uso delle serie temporali è rappresentato dalla semplicità nella rappresentazione grafica. Usando il dataset di esempio Nile presente in R, che è già un oggetto serie temporale – e posso verificarlo con il comando is.ts(Nile) -, mi basterà il comando:
plot.ts(Nile)
per avere un grafico dell’andamento nel tempo della mia variabile. Ovviamente posso usare i vari attributi xlab, ylab,main ecc… per rendere ancora più significativo e chiaro il mio grafico:
plot(Nile, xlab="Anno", ylab="Flusso annuale del fiume Nilo", main="Serie Storica di esempio Nilo")
Se il dataset contiene più di un oggetto serie temporale, potrò ottenere il grafico dei vari oggetti.
Usiamo un altro dataset di esempio presente in R chiamato EuStockMarkets che contiene, come si può facilmente immaginare, le quotazioni dei listini FTSE, CAC, DAX e SMI :
plot(EuStockMarkets, plot.type = "multiple",)
Oppure potrò disegnare i vari oggetti serie temporale tutti insieme – con colori differenti – nello stesso grafico:
plot.ts(EuStockMarkets, plot.type = "single", col=c("red","black","blue","green")) Se andiamo a rappresentare graficamente una serie temporale, noteremo quasi sempre tutta una serie di piccole variazioni che possono rendere assai arduo l’obiettivo di individuare tendenze importanti e fare previsioni per il futuro. Proprio per approcciare questo problema sono state sviluppate varie tecniche di «lisciamento» (smoothing), che possiamo per semplicità suddividere in due grandi famiglie: le tecniche che si basano su medie mobili e le tecniche esponenziali.
MEDIA MOBILE
Al posto del dato relativo al mese X calcolo la media di un numero n di mesi di cui X è il punto centrale.
La componente casuale, si compensa se mettiamo assieme diversi mesi, la sua media è uguale a 0 per un numero ragionevole di periodi.
La componente stagionale si ripete regolarmente nel corso dell’anno, allora se distribuisco l’effetto stagionale su tutti i 12 mesi, l’effetto scompare.
Con la media mobile ottengo tutte e due gli effetti voluti: compenso la casualità e “distribuisco” la stagionalità.
R, come vedremo a breve, ci fornisce un aiuto fondamentale mettendoci a disposizione tutta una serie di strumenti per effettuare le nostre analisi con la massima praticità.
Come abbiamo visto, possiamo creare facilmente una serie storica in R usando il comando di base ts().
Dal momento che ho intenzione di usare l’analisi delle serie storiche a fini SEO per ricavare una tendenza e fare una previsione, ho bisogno di importare i dati della mia vista Google Analytics in R.
Posso farlo in maniera “automatica”, usando l’utilissima libreria googleAnalyticsR (che magari sarà oggetto di un post successivo più dettagliato) oppure esportando i dati che mi servono via Query Explorer.
Vediamo il primo caso:
# uso la libreria googleAnalyticsR
# che ovviamente devo aver
# correttamente installato
library(googleAnalyticsR)
# setto l'id della mia vista
# codice numerico ids che trovo
# in Analytics o in Query Esplorer:
view_id <- xxxxxxxxxx
# Autorizzo Google Analytics
ga_auth()
# Recupero i dati che mi servono
# le sessioni da inizio 2017
# a tutto il 2019 - da Google Analytics
gadata <- google_analytics(view_id,
date_range = c("2017-01-01", "2019-12-31"),
metrics = "sessions",
dimensions = c("yearMonth"),
max = -1)
# Converto i dati in una serie temporale
# con cadenza mensile indicando frequency=12:
ga_ts <- ts(gadata$sessions, start = c(2017,01), end = c(2019,12), frequency = 12) La procedura per ottenere lo stesso risultato partendo da un file tsv è altrettanto semplice.
Per prima cosa vado sull’interfaccia di Query Explorer
e seleziono il mio account e la vista.
Scelgo la data di inizio e quella di fine (nel mio esempio dal primo gennaio 2017 al 31 dicembre 2019), la metrica che mi interessa (ga:sessions), la dimensione (ga:date).
Un clic sul bottone “Run Query” ed avrò a video il risultato, che posso scaricare in locale cliccando su “Download Results as TSV”.
Nel mio esempio salvo il file come nomefile.tsv
Ora apro con un editor di testo il file ed elimino le prime righe con le informazioni generali e in fondo al file il totale generale.
Se voglio, posso rinominare la riga di header, sostituendo “ga:date” con “data” e “ga:sessions” con “sessioni”, per una migliore leggibilità.
Non ci resta che importare il tsv e creare la serie temporale. Questione di due righe:
# Importo un semplicissimo dataset
# con data e utenti per giorno
sitomese <- read.csv('c:/percorso/nomefile.tsv', header = TRUE, sep = "\t")
# Converto i dati in una serie temporale
sitomese_ts <- ts(sitomese$sessioni, start = c(2017,01), end = c(2019,12), frequency = 12) Ora che ho la mia serie temporale, ho a disposizione una molteplicità di pacchetti R che mi forniscono tutti gli strumenti utili per analisi di qualunque tipo, dalle più elementari a quelle più approfondite.
Installo il pacchetto forecast per poter usare l’utilissima funzione ma():
library(forecast) sitomese.filt<-ma(sitomese_ts,order=12) sitomese.filt
In questo modo si è applicata una media ponderata alla nostra serie temporale, in modo da limitare l’effetto della stagionalità. Posso ora visualizzare la tendenza stimata con il sistema della media mobile:
lines(sitomese.filt,col="red")
Usando la differenza con la funzione diff() e un lag appropriato, ho la possibilità di eliminare il trend stagionale, qualora presente. Nel caso di una serie storica con dati mensili che comprende più anni di osservazioni, mi basta usare un lag=12, come in questo banale esempio nel quale ipotizzo di lavorare su una serie temporale x:
# elimino il trend stagionale # ipotizzo dati mensili e la presenza di stagionalità dx <- diff(x, lag=12) ts.plot(dx, main="Tendenza utenti destagionalizzata")
Abbiamo visto nel corso dell’articolo che il metodo classico di analisi delle serie storiche individua quattro influenze, o componenti e che lo scopo principale è proprio quello di scomporre la serie, per isolare le influenze delle varie componenti che determinano i valori della serie storica.
Passando dalla teoria alla pratica, vediamo come possiamo procedere.
Posso sfruttare la funzione decompose() del pacchetto stats per operare una classica decomposizione della mia serie temporale nelle componenti usando il sistema delle medie mobili e rappresentando in un solo chiarissimo grafico il tutto:
# decompongo la serie temporale. Ho scelto una decomposizione # moltiplicativa, ovviamente avrei # potuto scegliere una additiva componenti <- decompose(sitomese_ts, type ="multiplicative") names(componenti) # esploro nel grafico le componenti della serie storica plot(componenti)
Un’alternativa più raffinata per la decomposizione della serie è quella che usa il metodo loess (Locally Weighted Smoothing). Si tratta di un insieme di metodi non-parametrici che adattano i modelli di regressione polinomiale ai sottoinsiemi di dati. Utilizziamo a tale scopo la funzione stl() del pacchetto stats:
# uso stl per un decompose di tipo LOESS sitomese-loess <- stl(sitomese_ts, s.window="periodic") names(sitomese-loess) plot(sitomese-loess)
Le tecniche di lisciamento (smoothing) e previsione (forecast) ci offrono potenti modalità operative al fine di prevedere valori futuri dei dati delle serie temporali.
Al livello più elementare, il livellamento può essere realizzato utilizzando le medie mobili.
In R, possiamo usare HoltWinters, una funzione per eseguire il livellamento delle serie storiche.
La funzione contiene tre metodi di smoothing esponenziali. Tutti e tre i metodi usano la stessa funzione, HoltWinters. Tuttavia, possiamo invocarli separatamente in base ai valori dei parametri alpha, beta e gamma.
Il livellamento esponenziale Holt-Winters fornisce previsioni attendibili solo se non è presente autocorrelazione nei dati della serie temporale, cosa che può essere verificata, come vedremo tra breve in pratica, con la funzione acf e con un test Box–Pierce o Ljung–Box.
Dopo aver creato un modello di previsione, dobbiamo infatti valutarlo per capire se rappresenta correttamente i dati. In maniera analoga ad un modello di regressione, possiamo usare per questo scopo i residui. Se i residui seguono una distribuzione di tipo rumore bianco (white noise), allora la sequenza (o errore) dei residui è generata da un processo di tipo stocastico. E quindi il nostro modello ben rappresenta la serie temporale.
Vediamo un esempio. Supponiamo di avere una serie storica x :
# Usiamo la funzione forecast per una previsione: prossimi 6 periodi futuro.pre <- forecast(x.pre, h=6) # stampiamo a video un riepilogo summary(futuro.pre) # disegniamo il grafico plot(futuro.pre) # disegniamo il grafico dei residui per stimare l’autocorrelazione acf(futuro.pre$residuals,na.action = na.pass) # facciamo un test di autocorrelazione Box.test(futuro.pre$residuals)
L’autocorrelazione ci indica se i termini di una serie storica dipendono dal proprio passato.
Se consideriamo una serie temporale x di lunghezza n, l’autocorrelazione di lag 1 può essere stimata come la correlazione della coppia di osservazioni (x[t], x[t-1]).
R ci mette a disposizione un comodo comando: acf().
Usando:
acf(x, lag.max = 1, plot = FALSE)
sulla serie x viene automaticamente calcolata l’autocorrelazione di grado -1.
di default il comando acf(x) traccia un grafico, che mostra due righe orizzontali blu tratteggiate, che rappresentano l’intervallo di confidenza al 95%.
La stima dell’autocorrelazione è indicata dall’altezza delle barre verticali (ovviamente l’autocorrelazione al grado 0 è sempre 1).
L’intervallo di confidenza è usato per determinare la significatività statistica dell’autocorrelazione.
Mostro a mo’ d’esempio l’output della funzione acf() sulla serie storica Nile fornita da R:
Se il coefficiente di autocorrelazione diminuisce e cade rapidamente tra i bordi, ciò significa che i residui seguono una distribuzione di tipo rumore bianco. Non c’è evidente autocorrelazione.
Al contrario, se i coefficienti sono sempre sopra o sotto il limite, ciò significa che i residui sono autocorrelati.
Un test di autocorrelazione Ljung-Box è una forma particolare di test delle ipotesi, e fornisce un valore p come output, valore che ci consente di capire se rigettare l’ipotesi nulla o meno.
Applichiamo la funzione box.test sulla sequenza residua; troviamo il valore p. Se esso è superiore al valore di α non possiamo rifiutare l’ipotesi nulla. Cioè, i residui sono rumore bianco e ciò dimostra che il nostro modello “funziona bene” nella previsione del valore.
Vediamo il tutto in azione nel nostro esempio SEO relativo al traffico di un sito web, usando anche la libreria highcharter per una migliore visualizzazione dell’output:
# Per prima cosa carico le librerie che mi servono
library(googleAnalyticsR)
# per leggere i dati Analytics
library(forecast)
# per le previsioni su serie temporali
library(highcharter)
# per ottenere il grafico
# qua inserisco il codice ID
# per trovarlo basta loggarsi
# in Analytics Query Explorer
# https://ga-dev-tools.appspot.com/query-explorer/
# e leggere il valore "ids"
# per la vista che ci interessa
view_id <- xxxxxxxxx
# Autorizzo Google Analytics
ga_auth()
# e poi recupero i dati da Google Analytics
sitomese <- google_analytics_4(view_id,
date_range = c("2017-01-01", "2019-12-31"),
metrics = "sessions",
dimensions = c("yearMonth"),
max = -1)
# nb: la dimensione dei miei dati è annomese
# Ora esprimo i dati come serie temporale
sitomese_ts <- ts(sitomese$sessions, start = c(2017,01), end = c(2019,12), frequency = 12)
# Calcolo il livellamento Holt-Winters
previsione <- HoltWinters(sitomese_ts)
# Genero una previsione per i prossimi 12 mesi
hchart(forecast(previsione, h = 12)) A chi legge il compito di testare la bontà del modello previsionale.
L’uso del metodo di livellamento esponenziale richiede che i residui non siano correlati. Nei casi reali, questo è abbastanza improbabile. Abbiamo però altri strumenti a disposizione per affrontare anche questi casi: R ci fornisce la funzione ARIMA per costruire modelli di serie temporale che prendano in considerazione l’autocorrelazione.
L’utilissima funzione arima.sim consente di simulare un processo ARIMA generando i dati di una serie temporale ad hoc.
Attraverso questa funzione, dunque, possiamo iniziare a vedere due modelli di serie temporali basilari: il rumore bianco (white noise) e la passeggiata aleatoria (random walk).
Un modello ARIMA consta di tre componenti: ARIMA(p,d,q).
Il rumore bianco è l’esempio più elementare di processo stazionario. Le caratteristiche salienti:
Il modello white noise in termini ARIMA è dunque ARIMA(0,0,0).
Andiamo allora a simulare una serie temporale di questo tipo:
wn <- arima.sim(model = list(order = c(0,0,0)), n=100) ts.plot(wn)
Il Random Walk è un semplice esempio di processo non-stazionario. Presenta queste caratteristiche salienti:
Anche il modello della passeggiata aleatoria è un modello di base di serie temporale, e può essere facilmente simulato con la nostra funzione arima.sim.
Il modello Random Walk è la somma cumulativa di serie Rumore Bianco aventi media zero.
Da questo discende che la prima serie differenziata di una serie Random Walk è una serie White Noise!
Il modello ARIMA per una serie Random Walk è ARIMA(0,1,0).
Generiamo una serie di questo tipo e visualizziamola:
RW <- arima.sim(model = list(order = c(0,1,0)), n = 100) ts.plot(RW)
vediamo la controprova di quanto affermato sopra:
RWdiff <- diff(RW) ts.plot(RW-diff)
Otteniamo proprio una serie di tipo White Noise:
Il modello autoregressivo integrato a media mobile (ARIMA – AutoRegressive Integrated Moving Average) è noto anche come modello Box-Jenkins, dal nome degli statistici George Box e Gwilym Jenkins.
Lo scopo di ARIMA è di trovare il modello che meglio rappresenti i valori di una serie storica.
Un modello ARIMA può essere espresso come ARIMA (p, d, q), dove, come abbiamo già visto, p è l’ordine del modello autoregressivo, d indica il grado di differenziazione e q indica l’ordine della media mobile..
Operativamente, possiamo definire cinque passaggi per adattare le serie storiche a un modello ARIMA:
1. Simulo un processo ARIMA grazie alla funzione arima.sim() e disegno il grafico:
simts <- arima.sim(list(order = c(1,1,0), ar = 0.64), n = 100) plot(simts)
2. Faccio la differenza per ottenere una serie storica stazionaria e traccio il grafico:
simts.diff <- diff(simts) plot(simts.diff)
3. Uso la funzione auto.arima per stimare i migliori valori di p,d, e q:
auto.arima(simts, ic="bic") Series: simts ARIMA(1,1,0) Coefficients: ar1 0.6331 s.e. 0.0760 sigma^2 estimated as 0.9433: log likelihood=-138.73 AIC=281.46 AICc=281.58 BIC=286.67
4. Creo il modello ARIMA con i valori p,d e q indicati (nel nostro esempio 1,1,0):
fit <- Arima(simts, order=c(1,1,0)) summary(fit)
5. Basandoci ora sul nostro modello ARIMA possiamo passare alla previsione di valori futuri della serie e a tracciare il grafico:
fit.prev <- forecast(fit) summary(fit.prev) plot(fit.prev)
Le aree colorate mostrano gli intervalli di confidenza all’80% e al 95%.
Passiamo infine a valutare la bontà del nostro modello con un grafico acf():
Come per il modello di livellamento esponenziale, possiamo usare la funzione acf per calcolare i residui e creare il diagramma di autocorrelazione. Poiché il coefficiente di autocorrelazione diminuisce rapidamente, i residui sono rumore bianco.
Possiamo anche fare un test Box-Pierce:
Box.test(fit.prev$residuals) Box-Pierce test data: fit.prev$residuals X-squared = 0.020633, df = 1, p-value = 0.8858
e otteniamo un valore p che indica la non rigettabilità dell’ipotesi nulla.
I test statistici sono strumenti fondamentali per l’analisi dei dati e la presa di decisioni informate. Scegliere…
Gli Alberi Decisionali sono un tipo di algoritmo di apprendimento automatico che utilizza una struttura…
Immaginiamo di voler trovare il percorso più veloce per raggiungere una destinazione in auto. Si…
Nel 1847, il matematico francese Augustin-Louis Cauchy stava lavorando su calcoli astronomici, quando ideò un…
La simulazione Monte Carlo è un metodo utilizzato per quantificare il rischio associato a un…
Abbiamo visto che la distribuzione binomiale si basa sull’ipotesi di una popolazione infinita N, condizione che si…