12.1 Stagionalità complessa

Finora abbiamo considerato soprattutto modelli stagionali relativamente semplici, adatti a dati trimestrali e mensili. Tuttavia, le serie temporali a frequenza più alta mostrano spesso andamenti stagionali più complicati. Per esempio, i dati giornalieri possono avere un modello settimanale così come un modello annuale. I dati orari hanno di solito tre tipi di stagionalità: giornaliera, settimanale e annuale. Anche i dati settimanali possono essere difficili da prevedere, dato che non c’è un numero intero di settimane in un anno, per cui il modello annuale ha un periodo stagionale di \(365.25/7 \approx 52.179\) in media. La maggior parte dei metodi che abbiamo considerato finora non è in grado di gestire queste complessità stagionali.

Non vogliamo necessariamente includere tutti i possibili periodi stagionali nei nostri modelli, ma solo quelli che probabilmente saranno presenti nei dati. Per esempio, se abbiamo solo 180 giorni di dati, possiamo ignorare la stagionalità annuale. Se i dati sono misure di un fenomeno naturale (per esempio, la temperatura), possiamo probabilmente ignorare con tranquillità la stagionalità settimanale.

La figura 12.1 mostra il numero di chiamate a una banca commerciale nordamericana con un intervallo di 5 minuti tra le 7:00 e le 21:05 di ogni giorno della settimana per un periodo di 33 settimane. Il pannello inferiore mostra le prime quattro settimane della stessa serie temporale. C’è una marcata stagionalità giornaliero di periodo 169 (ci sono 169 intervalli di 5 minuti al giorno), e una più debole stagionalità settimanale di periodo \(169 \times 5 = 845\). (I volumi di chiamata del lunedì tendono ad essere più alti del resto della settimana). Se fosse stata disponibile una serie più lunga di dati, avremmo potuto osservare anche un andamento stagionale annuale.

bank_calls %>%
  fill_gaps() %>%
  autoplot(Calls) +
  labs(y = "Chiamate",
       title = "Volume di chiamate alla banca ogni cinque minuti")
Volume di chiamate nell'arco di cinque minuti gestite nei giorni lavorativi di una settimana tra le 7:00 e le 21:05 in una granda banca commerciale dell'America del Nord. Pannello superiore: dati dal 3 Marzo al 24 Ottobre 2003. Pannello inferiore: prime quattro settimane di dati.

Figura 12.1: Volume di chiamate nell’arco di cinque minuti gestite nei giorni lavorativi di una settimana tra le 7:00 e le 21:05 in una granda banca commerciale dell’America del Nord. Pannello superiore: dati dal 3 Marzo al 24 Ottobre 2003. Pannello inferiore: prime quattro settimane di dati.

A parte le periodicità stagionali multiple, questa serie presenta un ulteriore elemento di complessità, dato dai valori mancanti tra i periodi lavorativi.

STL con periodicità stagionali multiple

La funzione STL() è progettata per trattare stagionalità multiple. Restituirà più componenti stagionali, oltre ad una componente di trend e una componente residua. In questo caso, dobbiamo reindicizzare lo tsibble per evitare i valori mancanti, e poi fornire esplicitamente i periodi stagionali.

calls <- bank_calls %>%
  mutate(t = row_number()) %>%
  update_tsibble(index = t, regular = TRUE)
calls %>%
  model(
    STL(sqrt(Calls) ~ season(period = 169) +
                      season(period = 5*169),
        robust = TRUE)
  ) %>%
  components() %>%
  autoplot() + labs(x = "Osservazione")
Decomposizione STL con stagionalità multipla per i dati sul volume delle chiamate.

Figura 12.2: Decomposizione STL con stagionalità multipla per i dati sul volume delle chiamate.

Appaiono due modelli stagionali, uno per l’ora del giorno (il terzo pannello), e uno per l’ora della settimana (il quarto pannello). Per interpretare correttamente questo grafico, è importante notare le scale verticali. In questo caso, il trend e la stagionalità settimanale hanno barre più larghe (e quindi intervalli relativamente più stretti) rispetto alle altre componenti, perché c’è poco trend nei dati, e la stagionalità settimanale è debole.

La decomposizione può anche essere utilizzata per le previsioni, con ciascuna delle componenti stagionali previste utilizzando un metodo naïve stagionale, e i dati destagionalizzati previsti utilizzando ETS.

Il codice è leggermente più complicato del solito perché dobbiamo aggiungere di nuovo gli identificativi temporali andati persi quando abbiamo reindicizzato lo tsibble per gestire i periodi con osservazioni mancanti. La trasformazione della radice quadrata usata nella decomposizione STL ha assicurato che le previsioni restino positive.

# Forecasts from STL+ETS decomposition
my_dcmp_spec <- decomposition_model(
  STL(sqrt(Calls) ~ season(period = 169) +
                    season(period = 5*169),
      robust = TRUE),
  ETS(season_adjust ~ season("N"))
)
fc <- calls %>%
  model(my_dcmp_spec) %>%
  forecast(h = 5 * 169)

# Add correct time stamps to fable
fc_with_times <- bank_calls %>%
  new_data(n = 7 * 24 * 60 / 5) %>%
  mutate(time = format(DateTime, format = "%H:%M:%S")) %>%
  filter(
    time %in% format(bank_calls$DateTime, format = "%H:%M:%S"),
    wday(DateTime, week_start = 1) <= 5
  ) %>%
  mutate(t = row_number() + max(calls$t)) %>%
  left_join(fc, by = "t") %>%
  as_fable(response = "Calls", distribution = Calls)

# Grafico dei risultati con i dati delle ultime 3 settimane
fc_with_times %>%
  fill_gaps() %>%
  autoplot(bank_calls %>% tail(14 * 169) %>% fill_gaps()) +
  labs(y = "Chiamate",
       title = "Volume di chiamate ogni cinque minuti alla banca")
#> Warning: Removed 100 rows containing missing values (`geom_line()`).
Previsioni dei dati sul volume delle chiamate usando una decomposizione STL in cui la previsione delle componenti stagionali usa un metodo stagionale naïve, e la previsione dei dati destagionalizzati usa ETS.

Figura 12.3: Previsioni dei dati sul volume delle chiamate usando una decomposizione STL in cui la previsione delle componenti stagionali usa un metodo stagionale naïve, e la previsione dei dati destagionalizzati usa ETS.

Regressione armonica dinamica con periodicità stagionali multiple

Con stagionalità multipla, possiamo usare i termini di Fourier come abbiamo fatto nei capitoli precedenti (si vedano le Sezioni 7.4 e 10.5). Poiché ci sono più stagionalità, abbiamo bisogno di aggiungere termini di Fourier per ogni periodo stagionale. In questo caso, i periodi stagionali sono 169 e 845, quindi i termini di Fourier sono della forma \[ \sin\left(\frac{2\pi kt}{169}\right), \quad \cos\left(\frac{2\pi kt}{169}\right), \quad \sin\left(\frac{2\pi kt}{845}\right), \quad \text{and} \quad \cos\left(\frac{2\pi kt}{845}\right), \] per \(k=1,2,\dots\). Come al solito, tali termini possono essere generati dalla funzione fourier().

Adatteremo un modello di regressione armonica dinamica con una struttura di errore ARIMA. Il numero totale di termini di Fourier per ogni periodo stagionale potrebbe essere selezionato minimizzando l’AICc. Tuttavia, per periodi stagionali elevati, questo tende a sovrastimare il numero di termini necessari, quindi useremo una scelta più soggettiva, con 10 termini per la stagionalità giornaliera e 5 per quella settimanale. Ancora una volta, useremo la trasformazione radice quadrata per garantire che le previsioni e gli intervalli di previsione rimangano positivi. Impostiamo \(D=d=0\) per gestire la non stazionarietà attraverso i termini di regressione e \(P=Q=0\) per gestire la stagionalità attraverso i termini di regressione.

fit <- calls %>%
  model(
    dhr = ARIMA(sqrt(Calls) ~ PDQ(0, 0, 0) + pdq(d = 0) +
                  fourier(period = 169, K = 10) +
                  fourier(period = 5*169, K = 5)))

fc <- fit %>% forecast(h = 5 * 169)

# Add correct time stamps to fable
fc_with_times <- bank_calls %>%
  new_data(n = 7 * 24 * 60 / 5) %>%
  mutate(time = format(DateTime, format = "%H:%M:%S")) %>%
  filter(
    time %in% format(bank_calls$DateTime, format = "%H:%M:%S"),
    wday(DateTime, week_start = 1) <= 5
  ) %>%
  mutate(t = row_number() + max(calls$t)) %>%
  left_join(fc, by = "t") %>%
  as_fable(response = "Calls", distribution = Calls)
# Plot results with last 3 weeks of data
fc_with_times %>%
  fill_gaps() %>%
  autoplot(bank_calls %>% tail(14 * 169) %>% fill_gaps()) +
  labs(y = "Chiamate",
       title = "Volume di chiamate ogni cinque minuti alla banca")
Previsioni da una regressione armonica dinamica applicata ai dati sul volume di chiamate.

Figura 12.4: Previsioni da una regressione armonica dinamica applicata ai dati sul volume di chiamate.

Questo è un modello grande, contenente 33 parametri: 4 coefficienti ARMA, 20 coefficienti di Fourier per il periodo 169, e 8 coefficienti di Fourier per il periodo 845. Non tutti i termini di Fourier per il periodo 845 sono usati perché c’è qualche sovrapposizione con i termini del periodo 169 (poiché \(845 = 5 \times 169\)).

Esempio: La domanda di energia elettrica

Un’applicazione comune di tali modelli è la modellazione della domanda di energia elettrica. La figura 12.5 mostra la domanda semioraria di energia elettrica (MWh) nel Victoria, Australia, durante il 2012-2014, insieme alle temperature (gradi Celsius) per lo stesso periodo per Melbourne (la più grande città del Victoria).

vic_elec %>%
  pivot_longer(Demand:Temperature, names_to = "Series") %>%
  ggplot(aes(x = Time, y = value)) +
  geom_line() +
  facet_grid(rows = vars(Series), scales = "free_y") +
  labs(y = "")
Domanda semioraria di energia elettrica e temperature corrispondenti nel 2012--2014, Victoria, Australia.

Figura 12.5: Domanda semioraria di energia elettrica e temperature corrispondenti nel 2012–2014, Victoria, Australia.

Il grafico della domanda di energia elettrica al variare della temperatura (Figura 12.6) mostra che c’è una relazione non lineare tra le due variabili, con la domanda che aumenta per le basse temperature (a causa del riscaldamento) e aumenta per le alte temperature (a causa del raffreddamento).

elec <- vic_elec %>%
  mutate(
    DOW = wday(Date, label = TRUE),
    Working_Day = !Holiday & !(DOW %in% c("Sat", "Sun")),
    Cooling = pmax(Temperature, 18)
  )
elec %>%
  ggplot(aes(x=Temperature, y=Demand, col=Working_Day)) +
  geom_point(alpha = 0.6) +
  labs(x="Temperatura (gradi Celsius)", y="Domanda (MWh)")
Domanda semioraria di energia elettrica per il Victoria, rappresentata graficamente rispetto alle corrispondenti temperature di Melbourne, la maggiore delle città del Victoria.

Figura 12.6: Domanda semioraria di energia elettrica per il Victoria, rappresentata graficamente rispetto alle corrispondenti temperature di Melbourne, la maggiore delle città del Victoria.

Adatteremo un modello di regressione con una funzione lineare a tratti della temperatura (contenente un nodo a 18 gradi), e termini di regressione armonica per tenere conto della stagionalità giornaliera. Ancora una volta, impostiamo gli ordini dei termini di Fourier soggettivamente, mentre usiamo l’AICc per selezionare l’ordine degli errori ARIMA.

fit <- elec %>%
  model(
    ARIMA(Demand ~ PDQ(0, 0, 0) + pdq(d = 0) +
          Temperature + Cooling + Working_Day +
          fourier(period = "day", K = 10) +
          fourier(period = "week", K = 5) +
          fourier(period = "year", K = 3))
  )

Prevedere con tali modelli è difficile, perché abbiamo bisogno di valori futuri delle variabili predittive. I valori futuri dei termini di Fourier sono facili da calcolare, ma le temperature future sono, ovviamente, sconosciute. Se siamo solo interessati a fare previsioni fino a una settimana in avanti, potremmo usare le previsioni di temperatura ottenute da un modello meteorologico. In alternativa, potremmo usare la previsione degli scenari (Sezione ??) e inserire gli eventuali modelli di temperatura. Nell’esempio seguente, abbiamo usato una ripetizione degli ultimi due giorni di temperature per generare possibili valori futuri della domanda.

elec_newdata <- new_data(elec, 2*48) %>%
  mutate(
    Temperature = tail(elec$Temperature, 2 * 48),
    Date = lubridate::as_date(Time),
    DOW = wday(Date, label = TRUE),
    Working_Day = (Date != "2015-01-01") &
                   !(DOW %in% c("Sat", "Sun")),
    Cooling = pmax(Temperature, 18)
  )
fc <- fit %>%
  forecast(new_data = elec_newdata)

fc %>%
  autoplot(elec %>% tail(10 * 48)) +
  labs(title="Domanda semioraria di energia elettrica: Victoria",
       y = "Domanda (MWh)", x = "Tempo [30m]")
Previsioni da un modello di regressione armonica applicato ai dati della domanda semioraria di energia elettrica.

Figura 12.7: Previsioni da un modello di regressione armonica applicato ai dati della domanda semioraria di energia elettrica.

Anche se le previsioni a breve termine sembrano ragionevoli, questo è un modello grezzo per un processo complicato. I residui, tracciati nella figura 12.8, dimostrano che ci sono molte informazioni che non sono state catturate con questo modello.

fit %>% gg_tsresiduals()
Diagnostiche sui residui per il modello di regressione armonica.

Figura 12.8: Diagnostiche sui residui per il modello di regressione armonica.

Versioni più sofisticate di questo modello che forniscono previsioni molto migliori sono descritte in Hyndman & Fan (2010) e Fan & Hyndman (2012).

Bibliografia

Fan, S., & Hyndman, R. J. (2012). Short-term load forecasting based on a semi-parametric additive model. IEEE Transactions on Power Systems, 27(1), 134–141. [DOI]
Hyndman, R. J., & Fan, S. (2010). Density forecasting for long-term peak electricity demand. IEEE Transactions on Power Systems, 25(2), 1142–1153. [DOI]