9.7 Modellazione ARIMA con fable

Come lavora ARIMA()?

La funzione ARIMA() nel pacchetto fable usa una variazione dell’algoritmo di Hyndman-Khandakar (Hyndman & Khandakar, 2008), che combina test di radici unitarie, minimizzazione dei criteri di AICc e MLE per ottenere un modello ARIMA. Gli argomenti di ARIMA() forniscono molte variazioni dell’algoritmo. Quello descritto qui è il comportamento di default.

Algoritmo di Hyndman-Khandakar per la scelta automatica di un modello ARIMA
  1. Il numero di differenze \(0 \le d\le 2\) viene determinato usando test KPSS ripetuti.
  1. I valori di \(p\) e \(q\) sono scelti minimizzando il criterio AICc dopo aver differenziato i dati \(d\) volte. Anzichè considerare ogni possibile combinazione di \(p\) e \(q\), l’algoritmo usa una procedura a passi per semplificare la ricerca.
  1. Inizialmente vengono adattati i seguenti quattro modelli:
    • ARIMA\((0,d,0)\),
    • ARIMA\((2,d,2)\),
    • ARIMA\((1,d,0)\),
    • ARIMA\((0,d,1)\).
    La costante è inclusa a meno che \(d=2\). Se \(d \le 1\), si considera anche il modello aggiuntivo:
    • ARIMA\((0,d,0)\) senza costante.
  1. Il miglior modello (quello col più piccolo valore di AICc) identificato al passo (a) viene impostato come “modello corrente”.
  1. Si considerano, quindi, le seguenti variazioni al modello corrente:
    • vary \(p\) e/o \(q\) del modello corrente \(\pm1\);
    • inclusione/esclusione di \(c\) dal modello corrente.
    Il miglior modello fra quelli considerati finora (il modello corrente o una di queste varianti) diventa il nuovo modello corrente.
  1. Si ripete il passo 2(c) fino a che non è più possibile ottenere valori inferiori del criterio AICc.
Un esempio illustrativo della procedura di  Hyndman-Khandakar

Figura 9.11: Un esempio illustrativo della procedura di Hyndman-Khandakar

La Figura 9.11 illustra graficamente, tramite un esempio, come funziona l’algoritmo di Hyndman-Khandakar per determinare automaticamente gli ordini di un modello ARMA. La griglia considera tutte le combinazioni di ordini ARMA(\(p,q\)) partendo dall’angolo in alto a sinistra con un ARMA(\(0,0\)), con l’ordine AR che aumenta lungo l’asse verticale e l’ordine MA che aumenta lungo l’asse orizzontale.

Le celle arancioni mostrano l’insieme iniziale di modelli considerati dall’algoritmo. In questo esempio, il modello ARMA(2,2) ha il valore AICc più basso tra questi modelli. Questo modello viene chiamato “modello corrente” ed è individuato dal cerchio nero. L’algoritmo cerca quindi tra i modelli vicini, come mostrato dalle frecce blu. Se viene trovato un modello migliore, questo diventa il nuovo “modello corrente”. In questo esempio, il nuovo “modello corrente” è il modello ARMA(3,3). L’algoritmo continua in questo modo finché non si trova un modello migliore. In questo esempio il modello identificato è un modello ARMA(4,2).

La procedura automatica passa a un nuovo “modello corrente” non appena viene identificato un modello migliore, senza provare tutti i modelli vicini. La ricerca completa (fatta cioè su tutti i modelli) viene effettuata quando greedy=FALSE.

La procedura automatica utilizza alcune approssimazioni per accelerare la ricerca. Queste approssimazioni possono essere evitate con l’argomento approximation=FALSE. È possibile che il modello con AICc minimo non venga trovato a causa di queste approssimazioni, o a causa dell’uso della procedura stepwise. Se si usa l’argomento stepwise=FALSE si cercherà il modello migliore su un insieme molto più ampio di modelli. Per una descrizione completa degli argomenti si consulti il file di aiuto.

Procedura di modellazione

La seguente procedura fornisce un utile approccio generale, quando si vuole adattare un modello ARIMA ad una serie storica (non stagionale).

  1. Rappresentare graficamente la serie e identificare le osservazioni anomale.
  2. Se necessario, trasformare i dati (usando una trasformata di Box-Cox) per stabilizzare la varianza.
  3. Se i dati sono non stazionari, differenziare i dati un numero di volte sufficiente a rendere stazionaria la serie.
  4. Esaminare i grafici delle funzioni di autocorrelazione empiriche, ACF/PACF: può essere appropriato un modello ARIMA(\(p,d,0\)) o ARIMA(\(0,d,q\))?
  5. Provare il modello o i modelli scelti e usa l’AICc per cercare un modello migliore.
  6. Controllare i residui del modello scelto utilizzando l’ACF dei residui e facendo un test portmanteau dei residui. Se non sembrano white noise, provare un modello leggermente modificato.
  7. Una volta che i residui sembrano essere white noise, usare il modello per prevedere.

L’algoritmo di Hyndman-Khandakar si occupa solo dei passi 3–5. Quindi, anche se lo si usa, bisognerà comunque occuparsi degli altri passi in autonomia.

Il processo è riassunto nella Figura 9.12.

Processo generale per prevedere usando un modello ARIMA.

Figura 9.12: Processo generale per prevedere usando un modello ARIMA.

Esempio: Esportazioni della Repubblica Centrafricana

Applichiamo la procedura appena descritta alla serie delle esportazioni della Repubblica Centrafricana riportate in Figura 9.13.

global_economy %>%
  filter(Code == "CAF") %>%
  autoplot(Exports) +
  labs(title="Central African Republic exports",
       y="% of GDP")
Esportazioni della Repubblica Centrafricana come percentuale del GDP.

Figura 9.13: Esportazioni della Repubblica Centrafricana come percentuale del GDP.

  1. Dall’osservazione del grafico della serie, si evince una certa non stazionarietà dovuta ad un trend decrescente. Il miglioramento che si verifica nel 1994 è dovuto a un nuovo governo che ha rovesciato la giunta militare e, inizialmente, ha avuto un certo successo, prima che i disordini causassero un ulteriore declino economico.

  2. Non sembra esserci alcuna evidenza di cambiamenti nella varianza, per cui non è necessaria alcuna trasformazione di Box-Cox.

  3. Per eliminare la non stazionarietà, applichiamo ai dati una differenza prima. La serie differenziata è mostrata nella Figura 9.14.

    global_economy %>%
      filter(Code == "CAF") %>%
      gg_tsdisplay(difference(Exports), plot_type='partial')
    Grafico, ACF e PACF per la serie delle differenze prime delle esportazioni della Repubblica Centrafricana.

    Figura 9.14: Grafico, ACF e PACF per la serie delle differenze prime delle esportazioni della Repubblica Centrafricana.

    La serie ora sembra essere stazionaria.

  4. La PACF riportata in Figura 9.14 suggerisce un modello AR(2); pertanto come candidato iniziale si può prendere un modello. la funzione ACF suggerisce un modello MA(3); quimndi un candidato alternativo è il mopdello ARIMA(0,1,3).

  5. Abbiamo adattato sia un modello ARIMA(2,1,0) che un modello ARIMA(0,1,3) insieme a due procedure di selezione automatica del modello, una che usa la procedura a passi descritta a passi in precedenza, e una che ricerca il modello migliore fra molti modelli.

    caf_fit <- global_economy %>%
      filter(Code == "CAF") %>%
      model(arima210 = ARIMA(Exports ~ pdq(2,1,0)),
            arima013 = ARIMA(Exports ~ pdq(0,1,3)),
            stepwise = ARIMA(Exports),
            search = ARIMA(Exports, stepwise=FALSE))
    
    caf_fit %>% pivot_longer(!Country, names_to = "Model name",
                             values_to = "Orders")
    #> # A mable: 4 x 3
    #> # Key:     Country, Model name [4]
    #>   Country                  `Model name`         Orders
    #>   <fct>                    <chr>               <model>
    #> 1 Central African Republic arima210     <ARIMA(2,1,0)>
    #> 2 Central African Republic arima013     <ARIMA(0,1,3)>
    #> 3 Central African Republic stepwise     <ARIMA(2,1,2)>
    #> 4 Central African Republic search       <ARIMA(3,1,0)>
    glance(caf_fit) %>% arrange(AICc) %>% select(.model:BIC)
    #> # A tibble: 4 × 6
    #>   .model   sigma2 log_lik   AIC  AICc   BIC
    #>   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
    #> 1 search     6.52   -133.  274.  275.  282.
    #> 2 arima210   6.71   -134.  275.  275.  281.
    #> 3 arima013   6.54   -133.  274.  275.  282.
    #> 4 stepwise   6.42   -132.  274.  275.  284.

    I quattro modelli considerati hanno valori AICc quasi identici. Dei modelli identificati, il metodo che utilizza la ricerca completa ha trovato che un modello ARIMA(3,1,0) che ha il valore AICc più basso, seguito dai modelli ARIMA(2,1,0) e ARIMA(0,1,3) — gli ultimi due sono i modelli che abbiamo identificato dai grafici ACF e PACF. La procedura di selezione automatica stepwise, invece, ha identificato un modello ARIMA(2,1,2), che ha il più alto valore AICc fra i quattro modelli.

  6. La funzione di autocorrelazione, ACF, dei residui del modello ARIMA(3,1,0) mostra che tutte le autocorrelazioni sono dentro le bande di confidenza, indicando, quindi, che i residui si comportano come white noise.

    caf_fit %>%
      select(search) %>%
      gg_tsresiduals()
    Grafici dei residui per il modello ARIMA(3,1,0).

    Figura 9.15: Grafici dei residui per il modello ARIMA(3,1,0).

    Inoltre, un test portmanteau ritorna un valore elevato del p-value, confermando l’ipotesi che i residui siano white noise.

    augment(caf_fit) %>%
      filter(.model=='search') %>%
      features(.innov, ljung_box, lag = 10, dof = 3)
    #> # A tibble: 1 × 4
    #>   Country                  .model lb_stat lb_pvalue
    #>   <fct>                    <chr>    <dbl>     <dbl>
    #> 1 Central African Republic search    5.75     0.569
  7. Le previsioni ottenute dal modello selezionato sono riportate nella Figura 9.16.

    caf_fit %>%
      forecast(h=5) %>%
      filter(.model=='search') %>%
      autoplot(global_economy)
    Previsioni per la serie delle esportazioni della Repubblica Centrafricana.

    Figura 9.16: Previsioni per la serie delle esportazioni della Repubblica Centrafricana.

    Si noti che le previsioni ottenute sono molto simili a quelle che si otterrebbero con un random walk (equivalente a un ARIMA(0,1,0)). L’inclusione di ulteriori termini AR e MA non ha portato a differenze sostanziali in questo esempio, anche se gli intervalli di previsione sono molto più stretti che per un modello random walk.

Il termine costante in R

Un modello ARIMA non stagionale può essere scritto come \[\begin{equation} \tag{9.3} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d y_t = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] o equivalentemente come \[\begin{equation} \tag{9.4} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d (y_t - \mu t^d/d!) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] dove \(c = \mu(1-\phi_1 - \cdots - \phi_p )\) e \(\mu\) è la media di \((1-B)^d y_t\). Il pacchetto fable usa la parametrizzazione dell’equazione (9.3) mentre molte altre routine di R usano la rappresentazione fornita dall’equazione (9.4).

Pertanto, l’inclusione di un termine costante in un modello ARIMA non stazionario equivale a indurre un trend polinomiale di ordine \(d\) nelle previsioni. (Se la costante viene omessa, le previsioni includono un trend polinomiale di ordine \(d-1\)). Quando \(d=0\), abbiamo il caso speciale in cui \(\mu\) è la media di \(y_t\).

Per default, la funzione ARIMA() determina automaticamente se una costante deve essere inclusa. Per \(d=0\) o \(d=1\), la costante viene inclusa se migliora il valore AICc. Se \(d>1\) la costante viene sempre omessa, poiché una tendenza quadratica o di ordine superiore è particolarmente pericolosa in fase di previsione.

La costante può essere specificata includendo 0 o 1 nella formula del modello (come l’intercetta in lm()). Per esempio, per selezionare automaticamente un modello ARIMA con costante, si può usare ARIMA(y ~ 1 + ...). Allo stesso modo, una costante può essere esclusa usando ARIMA(y ~ 0 + ...).

La rappresentazione grafica delle radici caratteristiche

(Questa è una sezione avanzata e può essere evitata se desiderato.)

Possiamo riscrivere l’equazione (9.3) come \[\phi(B) (1-B)^d y_t = c + \theta(B) \varepsilon_t\] dove \(\phi(B)= (1-\phi_1B - \cdots - \phi_p B^p)\) è un polinomio di ordine \(p\) in \(B\) e \(\theta(B) = (1 + \theta_1 B + \cdots + \theta_q B^q)\) è un polinomio di ordine \(q\) in \(B\).

Le condizioni di stazionarietà per il modello sono che le \(p\) radici complesse di \(\phi(B)\) stiano al di fuori del cerchio di raggio unitario, e le condizioni di invertibilità che le \(q\) radici complesse di \(\theta(B)\) stiano fuori dal cerchio di raggio unitario. Quindi, per vedere se il modello soddisfa le condizioni di invertibilità e/o di stazionarietà, possiamo usare un grafico che riporta le radici in relazione al cerchio (complesso) di raggio unitario.

In realtà, è più facile considerare l’inverso delle radici, dal momento che queste dovrebbero trovarsi tutte all’interno del cerchio di raggio unitario. Questo si può fare facilmente in R. Per il modello ARIMA(3,1,0) adattato alle esportazioni della Repubblica Centrafricana, otteniamo la figura (fig:armaroots).

gg_arma(caf_fit %>% select(Country, search))
Radici caratteristiche inverse per il modello ARIMA(3,1,0) adattato alle esportazioni della Repubblica Centrafricana.

Figura 9.17: Radici caratteristiche inverse per il modello ARIMA(3,1,0) adattato alle esportazioni della Repubblica Centrafricana.

I tre punti arancioni nel grafico corrispondono alle radici del polinomio \(\phi(B)\). Si può osservare, che si trovano tutte all’interno del cerchio di raggio unitario unitario, come ci aspettavamo, perché fable assicura che il modello stimato sia stazionario e invertibile. Qualsiasi radice vicina al cerchio unitario può essere numericamente instabile, e il modello corrispondente potrebbe non essere buono per le previsioni.

La funzione ARIMA() non restituirà mai un modello con (l’inverso delle) radici che si trovano al di fuori del cerchio unitario. Anche i modelli selezionati automaticamente dalla funzione ARIMA() non conterranno radici vicine a uno.index{ARIMA@} Di conseguenza, può accadere di trovare modelli con un valore AICc migliore di quello restituito da ARIMA(), ma tali modelli saranno potenzialmente problematici.

Bibliografia

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(1), 1–22. [DOI]