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


/root/anaconda3/envs/covid_ok/lib/python3.8/site-packages/statsmodels/tsa/statespace/sarimax.py:963: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
  warn('Non-stationary starting autoregressive parameters'
/root/anaconda3/envs/covid_ok/lib/python3.8/site-packages/statsmodels/tsa/statespace/sarimax.py:975: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

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) su dati trasformati con Box-Cox.

In modello migliore è attualmente

$$ \mathrm{ARIMA}(2, 1, 2)(1, 0, 1)[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_{=2})(P_{=1},D_{=0},Q_{=1})[s_{=7}] $$

RICALCOLO IN CORSO...

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 non ottimale (inferiore a 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). Ciononostante i dati sono stati trasformati con metodo Box-Cox per evitare risultati inferiori a zero (vedi https://otexts.com/fpp2/limits.html e https://otexts.com/fpp2/transformations.html).

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.

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.

Per il modello SARIMA utilizzeremo pertanto:

  • la derivata prima $d=1$ dei nuovi positivi
  • una stagionalità di 7 giorni ($s=7$)

e procederemo alla scelta del miglior modello che minimizzi il Bayesian Information Criterion per grid-search con pmdarima.

Performing stepwise search to minimize bic
 ARIMA(2,1,2)(1,0,1)[7]             : BIC=543.329, Time=0.62 sec
 ARIMA(0,1,0)(0,0,0)[7]             : BIC=862.719, Time=0.01 sec
 ARIMA(1,1,0)(1,0,0)[7]             : BIC=655.531, Time=0.06 sec
 ARIMA(0,1,1)(0,0,1)[7]             : BIC=752.885, Time=0.06 sec
 ARIMA(2,1,2)(0,0,1)[7]             : BIC=inf, Time=0.38 sec
 ARIMA(2,1,2)(1,0,0)[7]             : BIC=641.046, Time=0.29 sec
 ARIMA(2,1,2)(2,0,1)[7]             : BIC=548.713, Time=1.12 sec
 ARIMA(2,1,2)(1,0,2)[7]             : BIC=548.749, Time=1.28 sec
 ARIMA(2,1,2)(0,0,0)[7]             : BIC=inf, Time=0.31 sec
 ARIMA(2,1,2)(0,0,2)[7]             : BIC=670.073, Time=0.79 sec
 ARIMA(2,1,2)(2,0,0)[7]             : BIC=600.613, Time=0.51 sec
 ARIMA(2,1,2)(2,0,2)[7]             : BIC=555.582, Time=1.42 sec
 ARIMA(1,1,2)(1,0,1)[7]             : BIC=544.178, Time=0.45 sec
 ARIMA(2,1,1)(1,0,1)[7]             : BIC=567.685, Time=0.37 sec
 ARIMA(3,1,2)(1,0,1)[7]             : BIC=549.489, Time=0.75 sec
 ARIMA(2,1,3)(1,0,1)[7]             : BIC=549.678, Time=0.74 sec
 ARIMA(1,1,1)(1,0,1)[7]             : BIC=565.721, Time=0.18 sec
 ARIMA(1,1,3)(1,0,1)[7]             : BIC=543.417, Time=0.51 sec
 ARIMA(3,1,1)(1,0,1)[7]             : BIC=559.304, Time=0.66 sec
 ARIMA(3,1,3)(1,0,1)[7]             : BIC=555.243, Time=0.79 sec
 ARIMA(2,1,2)(1,0,1)[7] intercept   : BIC=inf, Time=0.85 sec

Best model:  ARIMA(2,1,2)(1,0,1)[7]          
Total fit time: 12.199 seconds
{'maxiter': 50,
 'method': 'lbfgs',
 'order': (2, 1, 2),
 'out_of_sample_size': 0,
 'scoring': 'mse',
 'scoring_args': {},
 'seasonal_order': (1, 0, 1, 7),
 'start_params': None,
 'suppress_warnings': True,
 'trend': None,
 'with_intercept': False}

SARIMAX Results
Dep. Variable: nuovi_positivi No. Observations: 574
Model: SARIMAX(2, 1, 2)x(1, 0, [1], 7) Log Likelihood -247.055
Date: Sun, 19 Sep 2021 AIC 508.109
Time: 16:42:35 BIC 538.565
Sample: 02-24-2020 HQIC 519.989
- 09-19-2021
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 -0.4671 0.263 -1.774 0.076 -0.983 0.049
ar.L2 -0.2695 0.067 -4.028 0.000 -0.401 -0.138
ma.L1 0.0549 0.264 0.208 0.835 -0.463 0.572
ma.L2 0.0978 0.132 0.739 0.460 -0.162 0.357
ar.S.L7 0.9756 0.008 126.081 0.000 0.960 0.991
ma.S.L7 -0.6873 0.032 -21.573 0.000 -0.750 -0.625
sigma2 0.1362 0.005 24.931 0.000 0.125 0.147
Ljung-Box (Q): 84.12 Jarque-Bera (JB): 193.62
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.28 Skew: 0.16
Prob(H) (two-sided): 0.00 Kurtosis: 5.83


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
mu lo 99% up 99%
index
ar.L1 -0.467104 -1.206107 0.271898
ar.L2 -0.269515 -0.457347 -0.081683
ma.L1 0.054924 -0.686191 0.796039
ma.L2 0.097803 -0.273631 0.469237
ar.S.L7 0.975610 0.953889 0.997331
ma.S.L7 -0.687339 -0.776774 -0.597905

Il modello migliore è molto semplice e con buone caratteristiche.