SCIR Model¶

Step 1 : extracting data from JHU with covid-19 Python library¶

We will use two elements provided in the JHU daily timeseries, namely the confirmed cases and the recovered cases.

Confirmed cases are what we note $I(t)$ in the model
Recovered cases for JHU are daily recoveries. You need to integrate these data to get $R(t)$

Let's load some libraries :

import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import numpy as np
from covid.plotting import plot_country_cases, get_world_map, plot_country_pace
from covid.database import WorldCovidMatcher, CountryCollection

Let's connect to the database and create the matchers for confirmed and recovered cases (the last possible matcher would be deaths) :

cm_c = WorldCovidMatcher('confirmed')                                                                                                                                           
cm_c.build_database()
cm_r = WorldCovidMatcher('recovered')                                                                                                                                           
cm_r.build_database()

We pick-up a judicious country list after discussion with colleagues in I11 meeting room and extract $I$ and $R$ (note only $R$ is cumulative in the extraction)

country_list=['France','Italy','Belgium','Austria']
I = cm_c.get_cases(CountryCollection(country_list), normalized=True, cumulative=False)
R = cm_r.get_cases(CountryCollection(country_list), normalized=True, cumulative=True)

We will first try to use the following equation $$\frac{dR}{dt}=\alpha I$$ and plot $\frac{R}{\int Idt}$ to have an estimate of $\alpha$

fig = go.Figure()
II=0
for k in I.keys():
    II+=I[k]
    fig.add_trace(go.Scatter(x=R[k].index, y=I[k], mode="lines", legendgroup="group2",name=k))    

fig.add_trace(go.Scatter(x=R[k].index, y=II, mode="lines", legendgroup="group2",name='Sum'))    
fig.update_layout(title='Sum of I',
                   xaxis_title='Date',
                   yaxis_title='Sum of I')
fig.update_yaxes(type="log")
fig.show()

fig = go.Figure()
for k in I.keys():
    fig.add_trace(go.Scatter(x=R[k].index, y=R[k]/np.cumsum(I[k]), mode="lines", legendgroup="group2",
                                         name=k,))
fig.update_layout(title='Rough estimate of α',
                   xaxis_title='Date',
                   yaxis_title='α')
fig.show()

fig = go.Figure()
II=0
for k in I.keys():
    II+=I[k]
    fig.add_trace(go.Scatter(x=R[k].index, y=I[k], mode="lines", legendgroup="group2",name=k))    

fig.add_trace(go.Scatter(x=R[k].index, y=II, mode="lines", legendgroup="group2",name='Sum'))    
fig.update_layout(title='Sum of I',
                   xaxis_title='Date',
                   yaxis_title='Sum of I')
fig.update_yaxes(type="log")
fig.show()

Well... Not so bad for a very simple test.

Now let's figure out what we can find about $\beta$. We will try two options :

First, considering the fact that at short times $I/R=\beta/\alpha-1$. Note that in the exponential regime we also have $\int Idt/\int Rdt=\beta/\alpha-1$, but the data is less noisy.
Second, using the exponential trend valid as long as $S\simeq N$ (or $\tilde{S}\simeq1$ in non dimensional variables) : $I\simeq I_0e^{(\beta-\alpha)t}$

fig = go.Figure()
for k in I.keys():
    fig.add_trace(go.Scatter(x=R[k].index, y=I[k]/R[k], mode="lines", legendgroup="group2",
                                         name=k,))
fig.update_layout(title='Rough estimate of 𝛽/α-1 from I/R',
                   xaxis_title='Date',
                   yaxis_title='𝛽/α-1')
fig.show()

fig = go.Figure()
for k in I.keys():
    fig.add_trace(go.Scatter(x=R[k].index, y=np.cumsum(I[k])/np.cumsum(R[k]), mode="lines", legendgroup="group2",
                                         name=k,))
fig.update_layout(title='Rough estimate of 𝛽/α-1 from ∫𝐼𝑑t/∫R𝑑t',
                   xaxis_title='Date',
                   yaxis_title='𝛽/α-1')
fig.show()

fig = go.Figure()
for k in I.keys():
    fig.add_trace(go.Scatter(x=I[k].index,y=np.cumsum(I[k]), mode="lines", legendgroup="group2",
                                         name=k))
fig.add_trace(go.Scatter(x=[I[k].index[30],I[k].index[40]],y=[.1,.1*np.exp(10/2.7)], mode="lines", legendgroup="group2",
                                         name='𝛽-α=1/2.7'))    
fig.update_layout(title='Rough estimate of 𝛽-α',
                   xaxis_title='Date',
                   yaxis_title='∫𝐼𝑑t')
fig.update_yaxes(type="log")
fig.show()

%%javascript
//hack to fix export
require.config({
  paths: {
    d3: 'https://cdnjs.cloudflare.com/ajax/libs/d3/5.9.2/d3',
    jquery: 'https://code.jquery.com/jquery-3.4.1.min',
    plotly: 'https://cdn.plot.ly/plotly-latest.min'
  },

  shim: {
    plotly: {
      deps: ['d3', 'jquery'],
      exports: 'plotly'
    }
  }
});

1/2.7

0.37037037037037035