COVID-19 Data Analysis

Previsione a 14 giorni

Previsione a 14 giorni con modello SARIMA.

Max Pierini

NB: questo articolo è un progetto aperto e ancora in fase di revisione e completamento. Per proposte di collaborazione, scrivere a info@epidata.it


I nuovi casi giornalieri previsti a 14 giorni sono stimati utilizzando un modello di autoregressione SARIMA ovvero $\mathrm{ARIMA}(p,d,q)(P,D,Q)[s]$ sui nuovi positivi osservati, i cui parametri ottimali sono stati determinati per grid-search (vedi oltre per dettagli).

In modello migliore è attualmente

$$ \mathrm{ARIMA}(2, 1, 0)(1, 0, 2)[7] $$

I parametri potrebbero essere ulteriormente affinati con successive grid-search mirate.

Previsioni

Metodo

SARIMA

I modelli ARIMA vengono solitamente scritti in forma semplificata

$$ \mathrm{ARIMA}(p,d,q)(P,D,Q)[s] $$

dove

  • $(p,d,q)$ sono rispettivamente gli iperparametri di auto regressione (AR), derivata (I) e media mobile (MA) della porzione non stagionale

  • $(P,D,Q)[s]$ sono rispettivamente gli iperparametri di auto regressione (AR), derivata (I) e media mobile (MA) della porzione con stagionalità $s$

e che in forma di operatori polinomiali corrisponde a

$$ \phi_p(B)\Phi_P(B^s)(1-B)^d(1-B^s)^D x_t \\ = \\ \theta_q(B) \Theta_Q (B^s) \varepsilon_t $$

in cui

  • $x_t$ è l'osservazione della variabile $x$ al tempo $t$

  • $\varepsilon_t$ è il rumore bianco $\varepsilon$ al tempo $t$ ed è identicamente e indipendentemente distribuito normale con media 0 e varianza $\sigma^2$ ovvero $\varepsilon \overset{iid}{\sim} \mathcal{N}(0, \sigma^2)$

  • l'operatore di back-shift $B^i$ è tale per cui $B^i y_t = y_{t-i}$ e pertanto i termini $(1 - B^i)^j y_t$ sono la derivata di ordine $j$ di $y_{t-i}$ infatti dove ad esempio $j=1$ otteniamo $(1 - B^i) y_t = y_{t-i} - y_{t-i-1}$

  • gli operatori $\phi_p(\cdot)$ e $\theta_q(\cdot)$ sono rispettavemente operatori polinomiali di autoregressione AR e media mobile MA, di ordine $p$ e $q$, della porzione non stagionale

  • gli operatori $\Phi_P(\cdot)$ e $\Theta_Q(\cdot)$ sono rispettavemente operatori polinomiali di autoregressione AR e media mobile MA, di ordine $P$ e $Q$, della porzione stagionale

Un operatore polinomiale $O_i(B^j)$ corrisponde a

$$ O_i(B^j) = 1 - \sum_{k=1}^{i} O_k B^{kj} $$

dove $O_k$ è il coefficiente di $B^{kj}$.

Il modello da noi scelto (come vedremo più avanti) corrisponde a

$$ \mathrm{ARIMA}(2,1,0)(1,0,2)[7] $$

quindi

$$ \mathrm{ARIMA}(p_{=2},d_{=1},q_{=0})(P_{=1},D_{=0},Q_{=2})[s_{=7}] $$

Sostituendo $(p,d,q)(P,D,Q)[s]$ nell'equazione a operatori polinomiali otteniamo

$$ \phi_2(B)\Phi_1(B^7)(1-B)^1(1-B^7)^0 x_t \\ = \\ \theta_0(B) \Theta_2 (B^7) \varepsilon_t $$

che possiamo semplificare in

$$ \phi_2(B)\Phi_1(B^7)(1-B) x_t \\ = \\ \Theta_2 (B^7) \varepsilon_t $$

Svolgendo gli operatori polinomiali AR e MA delle porzioni non stagionale e stagionale otteniamo pertanto il processo

$$ \begin{align*} (1- \phi_1 B - \phi_2 B^2) \\ (1 - \Phi_1 B^7) \\ (1 - B) x_t = \\ (1 - \Theta_1B^7 - \Theta_2B^{14}) \varepsilon_t \end{align*} $$

ovvero

$$ \begin{align*} (- \Phi_1 \phi_2 B^{10} - \phi_1 \Phi_1 B^9 \\ + \phi_2 \Phi_1 B^9 + \phi_1 \Phi_1 B^8 \\ + \Phi_1 B^8 - \Phi_1 B^7 \\ + \phi_2 B^3 + \phi_1 B^2 \\ - \phi_2 B^2 - \phi_1 B - B + 1) x_t \\ = \\ (1 - \Theta_1B^7 - \Theta_2B^{14}) \varepsilon_t \end{align*} $$

Espandendo gli operatori di back-shift

$$ \begin{align*} (- \Phi_1 \phi_2 x_{t-10} - \phi_1 \Phi_1 x_{t-9} \\ + \phi_2 \Phi_1 x_{t-9} + \phi_1 \Phi_1 x_{t-8} \\ + \Phi_1 x_{t-8} - \Phi_1 x_{t-7} \\ + \phi_2 x_{t-3} + \phi_1 x_{t-2} \\ - \phi_2 x_{t-2} - \phi_1 x_{t-1} - x_{t-1} + 1) x_t \\ = \\ \varepsilon_t - \Theta_1 \varepsilon_{t-7} - \Theta_2 \varepsilon_{t-14} \end{align*} $$

da cui otteniamo l'equazione esplicita

$$ \begin{align*} x_t & = x_{t-1} \\ & + \phi_1 x_{t-1} - \phi_1 x_{t-2} \\ & + \phi_2 x_{t-2} - \phi_2 x_{t-3} \\ & + \Phi_1 x_{t-7} - \Phi_1 x_{t-8} \\ & - \Phi_1 \phi_1 x_{t-8} + \Phi_1 \phi_1 x_{t-9} \\ & - \Phi_1 \phi_2 x_{t-9} + \Phi_1 \phi_2 x_{t-10} \\ & + \varepsilon_t - \Theta_1 \varepsilon_{t-7} - \Theta_2 \varepsilon_{t-14} \end{align*} $$

che possiamo infine semplificare in

$$ \begin{align*} x_t & = x_{t-1} \\ & + \phi_1 ( x_{t-1} - x_{t-2} ) \\ & + \phi_2 ( x_{t-2} - x_{t-3} ) \\ & + \Phi_1 ( x_{t-7} - x_{t-8} ) \\ & - \Phi_1 \phi_1 ( x_{t-8} - x_{t-9} ) \\ & - \Phi_1 \phi_2 ( x_{t-9} - x_{t-10} )\\ & + \varepsilon_t - \Theta_1 \varepsilon_{t-7} - \Theta_2 \varepsilon_{t-14} \end{align*} $$

dove i rumori bianchi $\varepsilon_t$ sono identicamente e indipendentemente distribuiti normali con media 0 e varianza $\sigma^2$, ovvero $\varepsilon_t \overset{iid}{\sim} \mathcal{N}(0, \sigma^2)$.

Abbiamo pertanto solamente 6 coefficienti (vedi statsmodels).

I coefficienti della porzione non stagionale

  • $\phi_1$, coefficiente dell'autoregressione AR, corripondente ad ar.L1
  • $\phi_2$, coefficiente dell'autoregressione AR, corripondente ad ar.L2

I coefficienti della porzione stagionale

  • $\Phi_1$, coefficiente dell'autoregressione AR, corripondente ad ar.S.L7
  • $\Theta_1$, coefficiente della media mobile MA, corripondente ad ma.S.L7
  • $\Theta_1$, coefficiente della media mobile MA, corripondente ad ma.S.L14

La varianza del rumore bianco $\varepsilon_t$

  • $\sigma^2$ corripondente a sigma2

Notiamo tre caratteristiche dei dati di cui è necessario tener conto prima di procedere alla scelta del modello SARIMA:

  1. Eteroschedasticità: la varianza aumenta all'aumentare della media; il test di Breusch-Pagan infatti ritorna un $p$-value notevolmente inferiore al più alto livello di significatività (0.1)
  2. Non-stazionarietà: i dati non sono stazionari ma presentano un chiaro trend non lineare; il test di Dickey-Fuller infatti ritorna un $p$-value superiore al più alto livello di significatività (0.1)
  3. Stagionalità: è chiaramente visibile una variazione ciclica (stagionale) di 7 giorni

I modelli SARIMA possono gestire la non-stazionarietà utilizzando la derivata delle osservazioni (parametro $d$) e la stagionalità tramite i parametri $(P,D,Q)[s]$. Scegliamo di non risolvere l'eteroschedasticità ma di riternerla una caratteristica propria dei dati raccolti (vedi ResearchGate).

Verifichiamo che l'utilizzo della derivata prima delle osservazioni sia in grado di risolvere il trend non lineare.

Come notiamo, il test di Dickey-Fuller sulla derivata prima ritorna un $p$-value minore di 0.05 e anche il test di Breusch-Pagan per l'eteroschedasticità è ora più accettabile (<0.1), anche se non ottimale.

Confermiamo infine la stagionalità settimanale con un periodogramma dei nuovi positivi, da cui si nota la fondamentale (primo armonico) di periodo $7$ giorni, il secondo e terzo armonico di $3.5$ e $2.\bar{3}$ giorni.