#!/usr/bin/env python
# coding: utf-8

# # Correlazione temporale tra diversi dati Covid

# ### Abstract
# Nelle due ondate di Covid in Italia, è stato più volte evidenziato che "il numero dei decessi è l'ultimo a diminuire" rispetto a tutti gli indicatori legati alla pandemia. E lo stesso vale per l'andamento dei guariti e dimessi.
# 
# Vero, ma **di quanti giorni stiamo parlando** rispetto, ad esempio, ai dati sui **nuovi positivi**?
# 
# Rispondo a questa e altre domande con un'analisi delle **serie temporali** dei dati Covid.

# ### Caricamento e overview dati Covid (da repository Protezione Civile)

# In[1]:


import pandas as pd
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots


# In[2]:


repo_pc = 'https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/'
file_andam_naz = 'dpc-covid19-ita-andamento-nazionale.csv'
raw_data = pd.read_csv(f'{repo_pc}{file_andam_naz}')


# In[3]:


raw_data.columns


# Limito i dati alla sola *seconda ondata*, dal primo settembre 2020 in avanti

# In[4]:


raw_data = raw_data[raw_data['data']>='2020-09-01T17:00:00']
raw_data


# I dati giornalieri non sono direttamente disponibili (ad eccezione dei *nuovi positivi*).
# 
# Li calcolo puntualmente in relazione allo *shift* sul giorno precedente

# In[5]:


raw_data['ricoverati_con_sintomi_giorno'] = raw_data['ricoverati_con_sintomi'] - raw_data['ricoverati_con_sintomi'].shift(1)
raw_data['terapia_intensiva_giorno'] = raw_data['terapia_intensiva'] - raw_data['terapia_intensiva'].shift(1)
raw_data['deceduti_giorno'] = raw_data['deceduti'] - raw_data['deceduti'].shift(1)
raw_data['tamponi_giorno'] = raw_data['tamponi'] - raw_data['tamponi'].shift(1)
raw_data['casi_testati_giorno'] = raw_data['casi_testati'] - raw_data['casi_testati'].shift(1)
raw_data['dimessi_guariti_giorno'] = raw_data['dimessi_guariti'] - raw_data['dimessi_guariti'].shift(1)


# Controllo eventuali anomalie.

# In[6]:


fig = make_subplots(rows=4, cols=2)

fig.add_trace(
    go.Bar(name='ricoverati_con_sintomi_giorno', x=raw_data['data'], y=raw_data['ricoverati_con_sintomi_giorno']),
    row=1, col=1
)
fig.add_trace(
    go.Bar(name='terapia_intensiva_giorno', x=raw_data['data'], y=raw_data['terapia_intensiva_giorno']),
    row=1, col=2
)
fig.add_trace(
    go.Bar(name='nuovi_positivi_giorno', x=raw_data['data'], y=raw_data['nuovi_positivi']),
    row=2, col=1
)
fig.add_trace(
    go.Bar(name='tamponi_giorno', x=raw_data['data'], y=raw_data['tamponi_giorno']),
    row=2, col=2
)
fig.add_trace(
    go.Bar(name='casi_testati_giorno', x=raw_data['data'], y=raw_data['casi_testati_giorno']),
    row=3, col=1
)
fig.add_trace(
    go.Bar(name='deceduti_giorno', x=raw_data['data'], y=raw_data['deceduti_giorno']),
    row=3, col=2
)
fig.add_trace(
    go.Bar(name='dimessi_guariti_giorno', x=raw_data['data'], y=raw_data['dimessi_guariti_giorno']),
    row=4, col=1
)

fig.update_layout(height=600, width=800, title_text='Variazioni giornaliere dei principali indicatori')
fig.update_layout(legend=dict(
    orientation="h",
    yanchor="bottom",
    y=-.3,
    xanchor="center",
    x=.5
))
fig.show()


# I casi testati meritano uno zoom.

# In[7]:


fig = px.bar(raw_data, x='data', y='casi_testati_giorno')
fig.show()


# Quasi 900.000 casi testati il 5 dicembre sono sicuramente un recupero nei dati o un errore.
# Annullo questa misurazione.

# In[8]:


raw_data.at[raw_data.index[raw_data['data']=='2020-12-05T17:00:00'][0],'casi_testati_giorno'] = np.NaN


# ### Auto-correlazione: quali dati presentano una vistosa periodicità

# In[9]:


days = 25
new_cases = \
  pd.DataFrame({'giorni': range(days),\
                'autocor_tamponi_eseguiti': \
                [raw_data['tamponi_giorno'].corr(raw_data['tamponi_giorno'].shift(i)) \
                            for i in range(days)], \
                'autocor_casi_testati': \
                [raw_data['casi_testati_giorno'].corr(raw_data['casi_testati_giorno'].shift(i)) \
                            for i in range(days)], \
                'autocor_nuovi_positivi': \
                [raw_data['nuovi_positivi'].corr(raw_data['nuovi_positivi'].shift(i)) \
                            for i in range(days)], \
                'autocor_nuovi_ricoverati': \
                [raw_data['ricoverati_con_sintomi_giorno'].corr(raw_data['ricoverati_con_sintomi_giorno'].shift(i)) \
                            for i in range(days)],\
                'autocor_nuove_TI': \
                [raw_data['terapia_intensiva_giorno'].corr(raw_data['terapia_intensiva_giorno'].shift(i)) \
                            for i in range(days)],\
                'autocor_nuovi_decessi': \
                [raw_data['deceduti_giorno'].corr(raw_data['deceduti_giorno'].shift(i)) \
                            for i in range(days)]
               })
fig = px.line(new_cases, \
              x='giorni', \
              y=['autocor_casi_testati','autocor_tamponi_eseguiti',\
                 'autocor_nuovi_positivi','autocor_nuovi_ricoverati','autocor_nuove_TI','autocor_nuovi_decessi'],\
              line_shape = 'spline',\
              title=f'Auto-correlazione su tamponi eseguiti / casi testati / '\
              'nuovi positivi / ricoverati / TI / decessi')
fig.add_vline(x = 7, line_dash='dash', line_color='green', \
              annotation_text='una<br>settimana', annotation_position='bottom right')
fig.update_yaxes(title_text='<b>Auto-correlazione</b>')
fig.update_xaxes(title_text='<b>Giorni</b>')
fig.show()


# Si nota chiaramente che il numero di **tamponi eseguiti** e di **casi testati** ha una fortissima **correlazione settimanale** (picchi ogni 7/14/21/.. giorni). In maniera minore, ma lo stesso si verifica anche per i **nuovi positivi**.
# 
# I nuovi **ricoveri** e le nuove **terapie intensive non presentano periodicità**. Va detto che questi valori sono in realtà la *variazione* nel numero di attualmente ricoverati e di persone attualmente in TI (non è disponibile il numero dei nuovi ingressi in ospedale o nuovi ingressi in TI). 

# ### Cross-correlazione: qual è il "ritardo" tra decessi e nuovi positivi

# Studiando la cross-correlazione è possibile stimare il *ritardo* tra decessi e nuovi positivi.

# In[10]:


days = 50
cases_death_corr = \
  pd.DataFrame({'giorni': range(days),\
                'crosscor_decessi_nuovi_positivi': \
                    [raw_data['deceduti_giorno'].rolling(7,center=True).mean()\
                     .corr(raw_data['nuovi_positivi'].rolling(7,center=True).mean().shift(i))\
                     for i in range(days)], \
                'crosscor_decessi_nuovi_ricoverati': \
                    [raw_data['deceduti_giorno'].rolling(7,center=True).mean()\
                     .corr(raw_data['ricoverati_con_sintomi_giorno'].rolling(7,center=True).mean().shift(i))\
                     for i in range(days)], \
                'crosscor_decessi_nuove_TI': \
                    [raw_data['deceduti_giorno'].rolling(7,center=True).mean()\
                     .corr(raw_data['terapia_intensiva_giorno'].rolling(7,center=True).mean().shift(i))\
                     for i in range(days)]
               })
fig = px.line(cases_death_corr, \
              x='giorni', \
              y='crosscor_decessi_nuovi_positivi', \
              line_shape = 'spline', \
              title='Cross-correlazione tra decessi e nuovi positivi, in media mobile a 7 giorni')
fig.update_yaxes(title_text='<b>Cross-correlazione tra decessi e positivi</b>')
fig.update_xaxes(title_text='<b>Giorni</b>')
fig.add_vline(x = cases_death_corr['crosscor_decessi_nuovi_positivi'].idxmax(),\
              line_dash='dash', line_color='green', \
              annotation_text='massima<br>correlazione', annotation_position='bottom right')
fig.show()


# In[11]:


giorni_max_cor = cases_death_corr['crosscor_decessi_nuovi_positivi'].idxmax()
valore_max_cor = round(cases_death_corr['crosscor_decessi_nuovi_positivi'].max(),4)
print(f'La maggiore correlazione tra decessi e nuovi positivi è dopo {giorni_max_cor} giorni')
print(f'Ed è pari a {valore_max_cor} su un massimo teorico di 1 (correlazione perfetta)')


# In termini più intuitivi: **i numeri dei decessi di oggi sono legati fortissimamente ai positivi di 13 giorni fa**.
# 
# Sovrapponendo i due grafici (con uno shift di 13 giorni e uso di doppia scala sulle y), si vede molto bene questa correlazione.

# In[12]:


fig = make_subplots(specs=[[{'secondary_y': True}]])
fig.add_trace(
    go.Scatter(x=raw_data['data'], y=raw_data['nuovi_positivi'], \
               name='numero positivi del giorno', line_shape='spline'),
    secondary_y=False,
)
fig.add_trace(
    go.Scatter(x=raw_data['data'], y=raw_data['deceduti_giorno'].shift(-13), \
               name='numero decessi di 13 giorni dopo', line_shape='spline'),
    secondary_y=True,
)

fig.update_layout(
    title_text='Andamento positivi e decessi (shift di 13 giorni e uso di doppia scala sulle y)'
)

fig.update_xaxes(title_text='Data')
fig.update_yaxes(title_text='<b>Positivi</b>', range=[0, 50000], secondary_y=False)
fig.update_yaxes(title_text='<b>Decessi di 13 giorni dopo</b>', range=[0, 1000],secondary_y=True)

fig.show()


# ### Cross-correlazione: tempi di guarigione

# E' interessante anche studiare la cross-correlazione tra nuovi positivi e dimessi/guariti.

# In[13]:


days = 50
healing_corr = \
  pd.DataFrame({'giorni': range(days),\
                'healing_corr': \
                    [raw_data['dimessi_guariti_giorno'].rolling(7,center=True).mean().\
                     corr(raw_data['nuovi_positivi'].rolling(7,center=True).mean().shift(i))\
                     for i in range(days)]
               })
fig = px.line(healing_corr, \
              x='giorni', \
              y='healing_corr', \
              line_shape = 'spline', \
              title='Cross-correlazione tra dimessi e nuovi ricoverati con sintomi')
fig.update_yaxes(title_text='<b>Cross-correlazione tra dimessi e positivi</b>')
fig.update_xaxes(title_text='<b>Giorni</b>')
fig.add_vline(x = healing_corr['healing_corr'].idxmax(),\
              line_dash='dash', line_color='green', \
              annotation_text='massima<br>correlazione', annotation_position='bottom right')
fig.show()


# In[14]:


giorni_max_cor_healing = healing_corr['healing_corr'].idxmax()
valore_max_cor_healing = round(healing_corr['healing_corr'].max(),4)
print(f'La maggiore correlazione tra guariti e nuovi positivi è dopo {giorni_max_cor_healing} giorni')
print(f'Ed è pari a {valore_max_cor_healing} su un massimo teorico di 1 (correlazione perfetta)')


# In altri termini: **i numeri dei guariti e dimessi di oggi sono legati fortissimamente ai nuovi positivi di 21 giorni fa**.
# 
# Anche in questo caso un grafico shiftato e con doppio asse y mostra chiaramente il fenomeno.

# In[15]:


fig = make_subplots(specs=[[{'secondary_y': True}]])
fig.add_trace(
    go.Scatter(x=raw_data['data'], y=raw_data['nuovi_positivi'], \
               name='numero positivi del giorno', line_shape='spline'),
    secondary_y=False,
)
fig.add_trace(
    go.Scatter(x=raw_data['data'], y=raw_data['dimessi_guariti_giorno'].shift(-21), \
               name='numero dimessi e guariti/dimessi di 21 giorni dopo', line_shape='spline'),
    secondary_y=True,
)

fig.update_layout(
    title_text='Andamento positivi e guariti/dimessi (shift di 21 giorni e uso di doppia scala sulle y)'
)

fig.update_xaxes(title_text='Data')
fig.update_yaxes(title_text='<b>Positivi</b>', range=[0, 45000], secondary_y=False)
fig.update_yaxes(title_text='<b>Guariti/dimessi di 21 giorni dopo</b>', range=[0, 45000], secondary_y=True)

fig.show()