COVID-19 Data Analysis

Previsione a 14 giorni

Previsione a 14 giorni con modello SARIMA.

Max Pierini

Predicting the Spread of SARS-CoV-2 in Italian Regions: The Calabria Case Study, February 2020–March 2022 (https://doi.org/10.3390/diseases10030038)


/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.

/root/anaconda3/envs/covid_ok/lib/python3.8/site-packages/statsmodels/base/transform.py:101: RuntimeWarning: invalid value encountered in power
  y = np.power(lmbda * x + 1, 1. / lmbda)

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: 858
Model: SARIMAX(2, 1, 2)x(1, 0, [1], 7) Log Likelihood 484.215
Date: Thu, 30 Jun 2022 AIC -954.429
Time: 16:31:55 BIC -921.237
Sample: 02-24-2020 HQIC -941.713
- 06-30-2022
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 -0.5169 0.161 -3.220 0.001 -0.831 -0.202
ar.L2 -0.2158 0.041 -5.250 0.000 -0.296 -0.135
ma.L1 -0.0165 0.162 -0.101 0.919 -0.335 0.302
ma.L2 -0.0037 0.089 -0.042 0.967 -0.179 0.172
ar.S.L7 0.9632 0.012 83.006 0.000 0.940 0.986
ma.S.L7 -0.5687 0.023 -24.447 0.000 -0.614 -0.523
sigma2 0.0186 0.001 31.383 0.000 0.017 0.020
Ljung-Box (Q): 116.78 Jarque-Bera (JB): 832.86
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.63 Skew: 0.12
Prob(H) (two-sided): 0.00 Kurtosis: 7.85


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
mu lo 99% up 99%
index
ar.L1 -0.516906 -0.967453 -0.066359
ar.L2 -0.215770 -0.331129 -0.100412
ma.L1 -0.016477 -0.472343 0.439389
ma.L2 -0.003730 -0.254768 0.247307
ar.S.L7 0.963156 0.930585 0.995728
ma.S.L7 -0.568675 -0.633973 -0.503378

Il modello migliore è molto semplice e con buone caratteristiche.