Scattered Thoughts - Epidemic modeling

# This code wrangles the data from JHU

import pandas as pd
import numpy as np

import math

from scipy import signal

import plotly.graph_objects as go
import plotly.express as px

from scipy.stats import expon
from scipy.stats import gamma
from scipy.stats import weibull_min

from numpy.random import default_rng
rng = default_rng()
import plotly
import plotly.io as pio
from IPython.display import display, HTML

## Tomas Mazak's workaround
plotly.offline.init_notebook_mode()
display(HTML(
    #'<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script src="https://cdn.plot.ly/plotly-3.0.1.js" charset="utf-8"></script>'
))

pio.renderers.default = "plotly_mimetype+notebook_connected"

pio.templates.default = "plotly_dark"

import datetime


# Import confirmed cases
conf_df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')

#Import deaths data
deaths_df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv')

# Import recovery data
rec_df = pd.read_csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv')

#iso_alpha = pd.read_csv('https://raw.githubusercontent.com/jeffufpost/sars-cov-2-world-tracker/master/data/iso_alpha.csv', index_col=0, header=0).T.iloc[0]
iso_alpha = pd.read_csv('https://raw.githubusercontent.com/jeffufpost/sars-cov-2-world-tracker/refs/heads/master/project/app/data/iso_alpha.csv', index_col=0, header=0)

# Wrangle the data

#print("Wrangling data by country.......")
# Consolidate countries (ie. frenc dom tom are included in France, etc..)
conf_df = conf_df.groupby("Country/Region")
conf_df = conf_df.sum().reset_index()
conf_df = conf_df.set_index('Country/Region')

deaths_df = deaths_df.groupby("Country/Region")
deaths_df = deaths_df.sum().reset_index()
deaths_df = deaths_df.set_index('Country/Region')

rec_df = rec_df.groupby("Country/Region")
rec_df = rec_df.sum().reset_index()
rec_df = rec_df.set_index('Country/Region')

# Remove Lat and Long columns
conf_df = conf_df.iloc[:,3:]
deaths_df = deaths_df.iloc[:,3:]
rec_df = rec_df.iloc[:,3:]

# Convert country names to correct format for search with pycountry
conf_df = conf_df.rename(index={'Congo (Brazzaville)': 'Congo', 'Congo (Kinshasa)': 'Congo, the Democratic Republic of the', 'Burma': 'Myanmar', 'Korea, South': 'Korea, Republic of', 'Laos': "Lao People's Democratic Republic", 'Taiwan*': 'Taiwan', "West Bank and Gaza":"Palestine, State of"})
# Convert country names to correct format for search with pycountry
deaths_df = deaths_df.rename(index={'Congo (Brazzaville)': 'Congo', 'Congo (Kinshasa)': 'Congo, the Democratic Republic of the', 'Burma': 'Myanmar', 'Korea, South': 'Korea, Republic of', 'Laos': "Lao People's Democratic Republic", 'Taiwan*': 'Taiwan', "West Bank and Gaza":"Palestine, State of"})
# Convert country names to correct format for search with pycountry
rec_df = rec_df.rename(index={'Congo (Brazzaville)': 'Congo', 'Congo (Kinshasa)': 'Congo, the Democratic Republic of the', 'Burma': 'Myanmar', 'Korea, South': 'Korea, Republic of', 'Laos': "Lao People's Democratic Republic", 'Taiwan*': 'Taiwan', "West Bank and Gaza":"Palestine, State of"})

# Convert dates to datime format
conf_df.columns = pd.to_datetime(conf_df.columns, format='mixed')
deaths_df.columns = pd.to_datetime(deaths_df.columns, format='mixed')
rec_df.columns = pd.to_datetime(rec_df.columns, format='mixed')

# Create a per day dataframe
#print("Creating new per day dataframes......")
# Create per day dataframes for cases, deaths, and recoveries - by pd.DatafRame.diff
conf_df_pd = conf_df.diff(axis=1)
deaths_df_pd = deaths_df.diff(axis=1)
rec_df_pd = rec_df.diff(axis=1)

#print("Create infected dataframe = conf - deaths - recoveries")
inf_df = conf_df - deaths_df - rec_df

conf_df_pd.iloc[:,0] = 0
rec_df_pd.iloc[:,0] = 0
deaths_df_pd.iloc[:,0] = 0
inf_df.iloc[:,0] = 0

#print("Adding dataframes of 1st, 2nd, and 3rd derivatives of number of infected")
firstdev = inf_df.apply(np.gradient, axis=1)
seconddev = firstdev.apply(np.gradient)
thirddev = seconddev.apply(np.gradient)

#print("Create series of first date above 100 confirmed cases.....")
# Create a column containing date at which 100 confirmed cases were reached, NaN if not reached yet
fda100 = conf_df[conf_df > 100].apply(pd.Series.first_valid_index, axis=1)

# Create dataframe for probability plot
probevent = iso_alpha.join(inf_df)
probevent['prev'] = probevent.iloc[:,-1] / probevent['SP.POP.TOTL']

# This code is the lowpass filter 
def lowpass(x, fc=0.05, pad=0): # starts with cutoff freq at 0.05 but can be adjusted, should use 0.02 for France for example
  
  # Pad with last value
  if pad == 1:
    i=0
    while i < 100:
      x=np.append(x, (x[-1]+x[-2]+x[-3])/3)
      i=i+1

  fs = 1  # Sampling frequency
  t = np.arange(len(x.T)) #select number of days done in SEIR model
  signala = x.T

  #fc = 0.05  # Cut-off frequency of the filter
  w = fc / (fs / 2) # Normalize the frequency
  b, a = signal.butter(5, w, 'low')
  return signal.filtfilt(b, a, signala), w

# This code creates the impulse responses
days = np.arange(100)
cdf = pd.DataFrame({
    'T_Latent': gamma.cdf(days, 1.8,loc=0.9,scale=(5.2-1.8)/0.9), 
    'T_Infectious': weibull_min.cdf(days, 2.3,loc=2,scale=20.11)
    })
h_L = cdf.diff().T_Latent
h_I = cdf.diff().T_Infectious
h_L[0] = 0
h_I[0] = 0

# This code is for the iterative deconvolution
# Let's define an iteration function:
def iter_deconv(alpha, impulse_response, input_signal, delay, comparator):
  conv=signal.fftconvolve(impulse_response, input_signal, mode='full')
  correction=np.roll(comparator-conv[:len(comparator)], delay)
  input_signal=np.floor(lowpass(alpha*correction+input_signal, 0.05, 0)[0])
  input_signal[input_signal<0]=0
  return input_signal

# Define a function to return MSE between two signals as a measure of goodness of fit
def msecalc(A, B):
  return ((A - B)**2).mean(axis=0)

Motivation for write-up

This is the 8th part of a multi-part series blog post on modeling in epidemiology.

While we have studied in detail the stochastc SEIR model, we have not compared this to actual real-world data.

The goal of this 8th installment is to expand on the 7th blog post and use deconvolution techniques on real world data to estimate the various parameters.

Most notably, we want to calculate \(R\) from available data from JHU, let’s see how to go about this.

Data available

We will use the JHU data which is updated daily and displayed graphically on my tracker here.

A lot of countries have had difficulty reporting reliable data, but a few have done so rather well.

We will have a closer look at these contries:

Austria
France
Germany
Switzerland
United States

The JHU datasets include the following data:

Daily cumulative confirmed cases
Daily cumulative recoveries
Daily cumulative deaths

From which we can calculate the following:

Daily new confirmed cases: \(I_{pd}\)
Daily new recoveries
Daily new deaths
Current number of cases: \(I\)

Analysis

We have seen previously how we can get \(E_{pd}\) from deconvolution of \(I_{pd}\).
We also know that:

\[E_{pd} = \beta ~ I ~ \frac{S}{N}\] \[\leftrightarrow \beta = \frac{N ~ E_{pd}}{I ~ S}\]

And we know:

\[R_0 = \frac{\beta}{\gamma}\] and \[R=R_0\frac{S}{N}\]

\[E_{pd}=\gamma~R~I\] \[\leftrightarrow R = \frac{E_{pd}}{\gamma~I}\]

This means we can find \(R\) by deconvolution of \(I_{pd}\).

Our aim is to find this graphically in the real-world data.

Real-world analysis

Generic analysis

Let’s first have a look at the genral way we will go about the analysis (i.e. not country specific).

We have daily data \(I_{pd}\) and \(R_{pd}\)
We have our assumed \(T_L\) and \(T_I\)

Note, depending on the data we have, \(R_{pd}\) can be:

sum of deaths and recoveries
only deaths
only recoveries

Analysis steps:

The first thing we want to check is whether \(h_I\circledast I_{pd} [j]\) gives us something close to \(R_{pd}\) If not, why not? (what kind of \(R_{pd}\) should we use?)
Can we get an estimated \(E_{pd}\) by deconvolution of \(I_{pd}\) ? What cutoff frequency should we use ?
What can that tell us about \(R\) and \(\beta\) ?

Austria

Lets first have a look at data from Austria

Use fc = 0.05

1. Checking \(h_I\)

Code

fig = go.Figure(data=[    
    go.Bar(name='Ipd', x=conf_df_pd.loc['Austria'].index, y=conf_df_pd.loc['Austria']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=conf_df_pd.loc['Austria'].index, y=lowpass(conf_df_pd.loc['Austria'], 0.05, 1)[0]),
    go.Bar(name='Rpd', x=rec_df_pd.loc['Austria'].index, y=rec_df_pd.loc['Austria']),
    go.Scatter(name='Rpd=lowpass(Rpd)', x=rec_df_pd.loc['Austria'].index, y=lowpass(rec_df_pd.loc['Austria'], 0.05, 1)[0]),
    go.Scatter(name='Rpd=conv(Ipd)', x=conf_df_pd.loc['Austria'].index, y=signal.fftconvolve(h_I, conf_df_pd.loc['Austria'], mode='full'))
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Austria: Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see the actual \(R_{pd}\) leads \(h_I[j]\circledast I_{pd}[j]\) by about 5 days.

There is a 20 day lag between peak \(I_{pd}\) and \(h_I[j]\circledast I_{pd}[j]\) as is expected since \(E[T_I] = 20.11\).

But there is only a 15 day lag between peak \(I_{pd}\) and actual \(R_{pd}\).

There are a few possibilites for why this is the case, including that maybe we haven’t assumed the correct distribution for \(T_I\).

However there is another reason why this could be. Testing, especially in the early days, took time, and it took time for a patient showing symptoms before he could be tested. This may simply explain the 5 day difference.

2. Estimating \(E_{pd}\) by deconvolution of \(I_{pd}\)

Code

#Settting up for deconvolution of Ipd

#regularization parameter
alpha=2

# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(conf_df_pd.loc['Austria'], 0.05, 1)[0])
Ipd[Ipd<0]=0


# Pad with last value
i=0
while i < 100:
  Ipd=np.append(Ipd, Ipd[-1])
  i=i+1

# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()

# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess

# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 10000000)
mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Austria'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Austria'])]))

itercount=0
while mse[-1] < mse[-2]:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Austria'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Austria'])]))
  print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")

Iteration #1: MSE= 280097.036893302
Iteration #2: MSE= 146227.1215060887
Iteration #3: MSE= 98427.81506218898
Iteration #4: MSE= 73637.12813586244
Iteration #5: MSE= 59569.91412113941
Iteration #6: MSE= 50477.19976409924
Iteration #7: MSE= 44172.759377519884
Iteration #8: MSE= 39507.84001002058
Iteration #9: MSE= 35832.886697091635
Iteration #10: MSE= 32875.44754973414
Iteration #11: MSE= 30378.109945193886
Iteration #12: MSE= 28271.241453430164
Iteration #13: MSE= 26418.481502302486
Iteration #14: MSE= 24812.840293680038
Iteration #15: MSE= 23369.96252175532
Iteration #16: MSE= 22096.63295282359
Iteration #17: MSE= 20937.04493230752
Iteration #18: MSE= 19894.25464484861
Iteration #19: MSE= 18948.959639641937
Iteration #20: MSE= 18087.458410139705
Iteration #21: MSE= 17300.24920923025
Iteration #22: MSE= 16577.668006392993
Iteration #23: MSE= 15919.451464138658
Iteration #24: MSE= 15300.907044372325
Iteration #25: MSE= 14744.675161225703
Iteration #26: MSE= 14210.683741229928
Iteration #27: MSE= 13738.119685289568
Iteration #28: MSE= 13277.832201974481
Iteration #29: MSE= 12873.332911900452
Iteration #30: MSE= 12469.245964828155
Iteration #31: MSE= 12118.224107328258
Iteration #32: MSE= 11759.266569961683
Iteration #33: MSE= 11449.777699700899
Iteration #34: MSE= 11131.659957279466
Iteration #35: MSE= 10862.139479901662
Iteration #36: MSE= 10573.068671552284
Iteration #37: MSE= 10333.756345607035
Iteration #38: MSE= 10076.654598951065
Iteration #39: MSE= 9869.633801315156
Iteration #40: MSE= 9633.22719271893
Iteration #41: MSE= 9446.34059525184
Iteration #42: MSE= 9234.269853278713
Iteration #43: MSE= 9069.962204752354
Iteration #44: MSE= 8873.17786679817
Iteration #45: MSE= 8722.147391985161
Iteration #46: MSE= 8542.091515068189
Iteration #47: MSE= 8408.644590365482
Iteration #48: MSE= 8242.506167912048
Iteration #49: MSE= 8124.5120742979625
Iteration #50: MSE= 7972.987989366172
Iteration #51: MSE= 7860.598971287295
Iteration #52: MSE= 7719.458842114416
Iteration #53: MSE= 7621.323727533386
Iteration #54: MSE= 7490.83303769503
Iteration #55: MSE= 7397.583821951359
Iteration #56: MSE= 7275.316302386207
Iteration #57: MSE= 7189.32491143722
Iteration #58: MSE= 7075.608628838437
Iteration #59: MSE= 6997.154682468184
Iteration #60: MSE= 6891.669715931714
Iteration #61: MSE= 6820.000022095949
Iteration #62: MSE= 6721.931341894739
Iteration #63: MSE= 6656.7407297932405
Iteration #64: MSE= 6565.26035425735
Iteration #65: MSE= 6504.275624274783
Iteration #66: MSE= 6415.65849730165
Iteration #67: MSE= 6356.806568959921
Iteration #68: MSE= 6271.206123333204
Iteration #69: MSE= 6219.307514472439
Iteration #70: MSE= 6142.708483765
Iteration #71: MSE= 6091.546359145574
Iteration #72: MSE= 6016.050209112622
Iteration #73: MSE= 5970.812362169446
Iteration #74: MSE= 5899.285625677148
Iteration #75: MSE= 5856.977454190732
Iteration #76: MSE= 5789.3161213834255
Iteration #77: MSE= 5747.087950803278
Iteration #78: MSE= 5679.500338125397
Iteration #79: MSE= 5644.187795502721
Iteration #80: MSE= 5580.280680635593
Iteration #81: MSE= 5545.408603337125
Iteration #82: MSE= 5484.431089279458
Iteration #83: MSE= 5453.903477449285
Iteration #84: MSE= 5392.874348290648
Iteration #85: MSE= 5364.440269167762
Iteration #86: MSE= 5303.190820142333
Iteration #87: MSE= 5276.1559599142865
Iteration #88: MSE= 5213.2843928561315
Iteration #89: MSE= 5191.00685069678
Iteration #90: MSE= 5129.745895988969
Iteration #91: MSE= 5106.63148781604
Iteration #92: MSE= 5042.710735774585
Iteration #93: MSE= 5024.520050637905
Iteration #94: MSE= 4959.777967931582
Iteration #95: MSE= 4943.725600949425
Iteration #96: MSE= 4875.271873179566
Iteration #97: MSE= 4862.827566294938
Iteration #98: MSE= 4792.683670844245
Iteration #99: MSE= 4781.061053161295
Iteration #100: MSE= 4708.305556553567
Iteration #101: MSE= 4698.486640356369
Iteration #102: MSE= 4625.1221124404055
Iteration #103: MSE= 4618.192025470725
Iteration #104: MSE= 4543.737145194148
Iteration #105: MSE= 4537.3247396068155
Iteration #106: MSE= 4457.709296744921
Iteration #107: MSE= 4452.973657491515
Iteration #108: MSE= 4374.22356769199
Iteration #109: MSE= 4371.7699337178765
Iteration #110: MSE= 4290.38013436198
Iteration #111: MSE= 4288.1246576608555
Iteration #112: MSE= 4207.8772648029535
Iteration #113: MSE= 4207.635756498308
Iteration #114: MSE= 4126.245976598282
Iteration #115: MSE= 4128.613977751894 so we use the result of the previous iteration.

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Code

# We can keep going the iteration until lowest MSE

#change alpha if you like
#alpha=2

i=0
while i < 10:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  print(msecalc(Ipd[:len(conf_df_pd.loc['Austria'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Austria'])]))
  i=i+1

4045.719721082141
4050.0822986420108
3968.020949588893
3971.9486870943533
3892.0850995742844
3899.536417475003
3820.717378063294
3828.929375478409
3750.356287600738
3762.533035688165

Code

fig = go.Figure(data=[    
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Austria'].index, y=Enext),
    go.Scatter(name='Ipd=conv(deconv(Ipd))', x=inf_df.loc['Austria'].index, y=signal.fftconvolve(h_L, Enext, mode='full')),
    go.Bar(name='Ipd', x=inf_df.loc['Austria'].index, y=conf_df_pd.loc['Austria']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=inf_df.loc['Austria'].index, y=lowpass(conf_df_pd.loc['Austria'], 0.05, 1)[0])
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Austria: Actual } I_{pd} \text{ vs. convolution of deconvolution of } I_{pd}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see that our estimate for \(\hat{E}_{pd}\) must be close to the reality of \(E_{pd}\) as \(I_{pd}\) is almost identical to \(\hat{E}_{pd} \circledast h_L\).

Of course, this holds only as long as our estimate of \(h_L\) is close to reality.

3. \(\beta\) and \(R\) from \(E_{pd}\) and \(I\)

As described above:

\[R = \frac{E_{pd}}{\gamma~I}\]

Code

# Calculate R
gam = 1/15 # As we say gamma is 1/20.11
R = Enext[:len(inf_df.loc['Austria'])]*(1/gam)/inf_df.loc['Austria']

fig = go.Figure(data=[    
    go.Scatter(name='R', x=inf_df.loc['Austria'].index, y=R),
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Austria'].index, y=Enext),
    go.Scatter(name='Inf', x=inf_df.loc['Austria'].index, y=inf_df.loc['Austria']),
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'R',
    title={
        'text':r'$\text{Austria: R }$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

France

Better to use fc = 0.02 here

1. Checking \(h_I\)

Code

fig = go.Figure(data=[    
    go.Bar(name='Ipd', x=conf_df_pd.loc['France'].index, y=conf_df_pd.loc['France']),
    go.Scatter(name='Ipd=lowpass2(Ipd)', x=conf_df_pd.loc['France'].index, y=lowpass(conf_df_pd.loc['France'], 0.02, 1)[0]),
    go.Bar(name='Rpd', x=rec_df_pd.loc['France'].index, y=rec_df_pd.loc['France']),
    go.Scatter(name='Rpd=lowpass(Rpd)', x=rec_df_pd.loc['France'].index, y=lowpass(rec_df_pd.loc['France'], 0.05, 1)[0]),
    go.Scatter(name='Rpd=conv(Ipd)', x=conf_df_pd.loc['France'].index, y=signal.fftconvolve(h_I, conf_df_pd.loc['France'], mode='full'))
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{France: Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Data is not very good, reason for using fc = 0.02 and it helps below:

2. Estimating \(E_{pd}\) by deconvolution of \(I_{pd}\)

Code

#Settting up for deconvolution of Ipd

#regularization parameter
alpha=2

# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(conf_df_pd.loc['France'], 0.02, 1)[0])
Ipd[Ipd<0]=0


# Pad with last value
#i=0
#while i < 100:
#  Ipd=np.append(Ipd, Ipd[-1])
#  i=i+1

  # Pad with last value
i=0
while i < 100:
  Ipd=np.append(Ipd, (Ipd[-1]+Ipd[-2]+Ipd[-3])/3)
  i=i+1

# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()

# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess

# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 10000000)
mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['France'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['France'])]))

itercount=0
while mse[-1] < mse[-2]:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['France'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['France'])]))
  print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Iteration #1: MSE= 5117960.805173197
Iteration #2: MSE= 3499161.0414016345
Iteration #3: MSE= 2468765.918751717
Iteration #4: MSE= 1786386.8815557766
Iteration #5: MSE= 1323470.1559443003
Iteration #6: MSE= 1001401.8255418137
Iteration #7: MSE= 772354.3244014457
Iteration #8: MSE= 607010.4444506937
Iteration #9: MSE= 484862.3581576968
Iteration #10: MSE= 394188.8154134461
Iteration #11: MSE= 325000.97457803506
Iteration #12: MSE= 272538.38354350056
Iteration #13: MSE= 231236.16894034666
Iteration #14: MSE= 199391.9603553421
Iteration #15: MSE= 173542.9406322499
Iteration #16: MSE= 153351.63805547243
Iteration #17: MSE= 136400.27538110517
Iteration #18: MSE= 123014.6376307402
Iteration #19: MSE= 111405.96650838591
Iteration #20: MSE= 102164.86985448647
Iteration #21: MSE= 93883.8951162736
Iteration #22: MSE= 87250.13892510024
Iteration #23: MSE= 81093.55956665453
Iteration #24: MSE= 76158.05529130949
Iteration #25: MSE= 71435.47195192643
Iteration #26: MSE= 67647.28730219323
Iteration #27: MSE= 63906.12267511958
Iteration #28: MSE= 60915.54720078659
Iteration #29: MSE= 57885.05294406694
Iteration #30: MSE= 55464.40767825074
Iteration #31: MSE= 52955.59887159959
Iteration #32: MSE= 50963.76269131391
Iteration #33: MSE= 48849.74613627011
Iteration #34: MSE= 47183.84792588964
Iteration #35: MSE= 45382.09526294367
Iteration #36: MSE= 43979.24164634863
Iteration #37: MSE= 42433.27883460782
Iteration #38: MSE= 41229.26384994886
Iteration #39: MSE= 39882.45044731062
Iteration #40: MSE= 38849.55664430575
Iteration #41: MSE= 37667.64431271743
Iteration #42: MSE= 36767.98060766154
Iteration #43: MSE= 35731.45868521261
Iteration #44: MSE= 34946.7550807061
Iteration #45: MSE= 34029.742900730496
Iteration #46: MSE= 33342.894475716756
Iteration #47: MSE= 32522.862385817636
Iteration #48: MSE= 31916.45201976579
Iteration #49: MSE= 31187.81327085697
Iteration #50: MSE= 30647.55167140667
Iteration #51: MSE= 29995.472893250575
Iteration #52: MSE= 29517.392032814387
Iteration #53: MSE= 28925.45668634593
Iteration #54: MSE= 28497.159636377663
Iteration #55: MSE= 27965.963581730128
Iteration #56: MSE= 27584.20587020581
Iteration #57: MSE= 27100.03407350276
Iteration #58: MSE= 26757.618630965688
Iteration #59: MSE= 26321.781601950872
Iteration #60: MSE= 26009.520783036874
Iteration #61: MSE= 25613.659472379128
Iteration #62: MSE= 25335.91232337881
Iteration #63: MSE= 24971.66701703756
Iteration #64: MSE= 24719.682474591726
Iteration #65: MSE= 24387.491275182463
Iteration #66: MSE= 24155.452973964537
Iteration #67: MSE= 23849.91316848782
Iteration #68: MSE= 23638.102536841725
Iteration #69: MSE= 23356.283402406825
Iteration #70: MSE= 23161.54110279197
Iteration #71: MSE= 22901.012327520115
Iteration #72: MSE= 22723.265209073295
Iteration #73: MSE= 22484.071668611046
Iteration #74: MSE= 22318.13541324509
Iteration #75: MSE= 22097.284420563035
Iteration #76: MSE= 21945.821745274046
Iteration #77: MSE= 21739.42407249935
Iteration #78: MSE= 21599.915034412476
Iteration #79: MSE= 21408.719014736234
Iteration #80: MSE= 21278.935432505943
Iteration #81: MSE= 21105.77383027484
Iteration #82: MSE= 20986.225901249083
Iteration #83: MSE= 20822.080589428246
Iteration #84: MSE= 20710.126618486927
Iteration #85: MSE= 20553.066562937816
Iteration #86: MSE= 20449.062399898816
Iteration #87: MSE= 20305.357670963087
Iteration #88: MSE= 20209.30971064579
Iteration #89: MSE= 20073.969643279695
Iteration #90: MSE= 19982.680108911507
Iteration #91: MSE= 19859.311006546966
Iteration #92: MSE= 19772.48969804378
Iteration #93: MSE= 19649.49714898565
Iteration #94: MSE= 19570.370507858686
Iteration #95: MSE= 19458.989766479433
Iteration #96: MSE= 19383.155540388165
Iteration #97: MSE= 19276.965354173233
Iteration #98: MSE= 19205.128971845577
Iteration #99: MSE= 19105.679980529276
Iteration #100: MSE= 19037.709165334334
Iteration #101: MSE= 18943.756611711768
Iteration #102: MSE= 18881.742290247254
Iteration #103: MSE= 18792.71228377153
Iteration #104: MSE= 18733.09712431994
Iteration #105: MSE= 18649.356922243474
Iteration #106: MSE= 18594.00893076957
Iteration #107: MSE= 18512.684919970667
Iteration #108: MSE= 18459.506144369676
Iteration #109: MSE= 18385.092400136822
Iteration #110: MSE= 18335.618303809842
Iteration #111: MSE= 18263.364853763876
Iteration #112: MSE= 18215.161051658226
Iteration #113: MSE= 18146.82167886694
Iteration #114: MSE= 18101.289189910443
Iteration #115: MSE= 18035.783683764203
Iteration #116: MSE= 17989.563157832228
Iteration #117: MSE= 17928.81812671534
Iteration #118: MSE= 17886.32849765552
Iteration #119: MSE= 17826.681771837786
Iteration #120: MSE= 17786.231816142128
Iteration #121: MSE= 17730.421890696078
Iteration #122: MSE= 17694.008622428857
Iteration #123: MSE= 17636.49634457983
Iteration #124: MSE= 17600.93097532715
Iteration #125: MSE= 17547.703385247645
Iteration #126: MSE= 17511.03327740333
Iteration #127: MSE= 17460.993258333427
Iteration #128: MSE= 17426.289098588026
Iteration #129: MSE= 17377.923688590105
Iteration #130: MSE= 17346.45974787351
Iteration #131: MSE= 17300.301501174756
Iteration #132: MSE= 17271.053601704396
Iteration #133: MSE= 17226.252685626314
Iteration #134: MSE= 17196.263116138143
Iteration #135: MSE= 17154.036524512954
Iteration #136: MSE= 17126.639467699053
Iteration #137: MSE= 17086.975999132606
Iteration #138: MSE= 17061.369020511705
Iteration #139: MSE= 17020.043412831164
Iteration #140: MSE= 16993.086016782392
Iteration #141: MSE= 16956.122216952095
Iteration #142: MSE= 16932.041833439707
Iteration #143: MSE= 16894.63026082054
Iteration #144: MSE= 16871.222439960944
Iteration #145: MSE= 16833.920715595737
Iteration #146: MSE= 16810.430259674067
Iteration #147: MSE= 16774.382336436876
Iteration #148: MSE= 16752.297663731042
Iteration #149: MSE= 16718.592035531652
Iteration #150: MSE= 16694.73120244851
Iteration #151: MSE= 16662.216563360376
Iteration #152: MSE= 16642.24930920911
Iteration #153: MSE= 16609.24851833289
Iteration #154: MSE= 16589.42491903914
Iteration #155: MSE= 16557.844463090674
Iteration #156: MSE= 16539.764780656278
Iteration #157: MSE= 16508.916637074104
Iteration #158: MSE= 16491.7120233013
Iteration #159: MSE= 16462.236392887324
Iteration #160: MSE= 16445.840009513857
Iteration #161: MSE= 16417.842786595476
Iteration #162: MSE= 16401.060581952508
Iteration #163: MSE= 16372.176192160174
Iteration #164: MSE= 16355.619323797577
Iteration #165: MSE= 16329.053873744771
Iteration #166: MSE= 16312.78745927229
Iteration #167: MSE= 16287.525501496004
Iteration #168: MSE= 16273.514281171914
Iteration #169: MSE= 16248.574012362944
Iteration #170: MSE= 16233.502381220185
Iteration #171: MSE= 16207.6173230868
Iteration #172: MSE= 16193.646902063929
Iteration #173: MSE= 16168.933010768356
Iteration #174: MSE= 16156.30063039807
Iteration #175: MSE= 16131.72764167814
Iteration #176: MSE= 16119.023096060057
Iteration #177: MSE= 16095.666422294615
Iteration #178: MSE= 16082.610177954266
Iteration #179: MSE= 16058.69008359635
Iteration #180: MSE= 16046.820857961895
Iteration #181: MSE= 16022.796541476928
Iteration #182: MSE= 16011.06818158278
Iteration #183: MSE= 15987.791123471306
Iteration #184: MSE= 15975.712302600152
Iteration #185: MSE= 15953.552933811581
Iteration #186: MSE= 15943.269293004003
Iteration #187: MSE= 15920.136222033112
Iteration #188: MSE= 15909.691648290083
Iteration #189: MSE= 15889.131431687952
Iteration #190: MSE= 15878.226783985763
Iteration #191: MSE= 15855.829456398604
Iteration #192: MSE= 15847.398353308003
Iteration #193: MSE= 15826.666109702275
Iteration #194: MSE= 15817.655034636287
Iteration #195: MSE= 15796.53948321836
Iteration #196: MSE= 15788.543369824518
Iteration #197: MSE= 15768.170596623606
Iteration #198: MSE= 15759.75542018989
Iteration #199: MSE= 15739.275448574028
Iteration #200: MSE= 15731.227727365633
Iteration #201: MSE= 15711.730093467388
Iteration #202: MSE= 15704.817237528549
Iteration #203: MSE= 15685.762132986054
Iteration #204: MSE= 15678.633662420483
Iteration #205: MSE= 15660.579766904706
Iteration #206: MSE= 15653.292821935827
Iteration #207: MSE= 15636.430033246625
Iteration #208: MSE= 15627.554356135823
Iteration #209: MSE= 15610.945604029603
Iteration #210: MSE= 15603.746862946895
Iteration #211: MSE= 15586.156432399022
Iteration #212: MSE= 15579.773727488786
Iteration #213: MSE= 15561.478054020394
Iteration #214: MSE= 15555.714926266843
Iteration #215: MSE= 15538.27924638053
Iteration #216: MSE= 15531.54785418079
Iteration #217: MSE= 15514.882884944716
Iteration #218: MSE= 15508.316905943075
Iteration #219: MSE= 15492.249877602435
Iteration #220: MSE= 15485.058416271666
Iteration #221: MSE= 15468.206964663017
Iteration #222: MSE= 15462.830324317474
Iteration #223: MSE= 15446.324515331247
Iteration #224: MSE= 15441.452382089226
Iteration #225: MSE= 15423.787783538908
Iteration #226: MSE= 15417.845641950138
Iteration #227: MSE= 15401.215287285762
Iteration #228: MSE= 15396.874555379341
Iteration #229: MSE= 15379.768939585469
Iteration #230: MSE= 15374.703793746698
Iteration #231: MSE= 15359.754242944002
Iteration #232: MSE= 15354.447392813581
Iteration #233: MSE= 15338.743428470243
Iteration #234: MSE= 15333.713841523362
Iteration #235: MSE= 15317.894920520106
Iteration #236: MSE= 15312.24483244489
Iteration #237: MSE= 15296.98486764927
Iteration #238: MSE= 15292.349722769733
Iteration #239: MSE= 15277.906982419268
Iteration #240: MSE= 15273.720286443047
Iteration #241: MSE= 15258.711851375814
Iteration #242: MSE= 15252.20364513739
Iteration #243: MSE= 15238.402464047384
Iteration #244: MSE= 15233.25748777765
Iteration #245: MSE= 15219.778835759176
Iteration #246: MSE= 15216.151753967033
Iteration #247: MSE= 15202.019079582995
Iteration #248: MSE= 15197.347316535446
Iteration #249: MSE= 15184.84686249038
Iteration #250: MSE= 15180.47417161207
Iteration #251: MSE= 15167.559836384089
Iteration #252: MSE= 15162.56194815394
Iteration #253: MSE= 15148.816471879534
Iteration #254: MSE= 15144.256996856007
Iteration #255: MSE= 15131.150515850835
Iteration #256: MSE= 15126.616699819802
Iteration #257: MSE= 15113.681035811553
Iteration #258: MSE= 15108.980866302414
Iteration #259: MSE= 15096.749746363821
Iteration #260: MSE= 15092.951356259222
Iteration #261: MSE= 15081.483567077576
Iteration #262: MSE= 15077.14782112155
Iteration #263: MSE= 15064.31830976379
Iteration #264: MSE= 15059.95924016886
Iteration #265: MSE= 15048.758065753747
Iteration #266: MSE= 15045.17066159609
Iteration #267: MSE= 15032.866484121587
Iteration #268: MSE= 15030.628370265009
Iteration #269: MSE= 15018.632308669052
Iteration #270: MSE= 15015.045680301353
Iteration #271: MSE= 15004.19357507284
Iteration #272: MSE= 15000.58468076094
Iteration #273: MSE= 14987.897025856882
Iteration #274: MSE= 14985.463060832295
Iteration #275: MSE= 14974.554599073244
Iteration #276: MSE= 14971.474782842488
Iteration #277: MSE= 14960.29354617457
Iteration #278: MSE= 14957.14625752821
Iteration #279: MSE= 14946.118517691404
Iteration #280: MSE= 14944.113964712637
Iteration #281: MSE= 14933.06954366123
Iteration #282: MSE= 14929.630066180982
Iteration #283: MSE= 14919.555775191899
Iteration #284: MSE= 14916.466318904275
Iteration #285: MSE= 14906.397882291029
Iteration #286: MSE= 14903.295906239666
Iteration #287: MSE= 14892.426019899116
Iteration #288: MSE= 14889.512614366968
Iteration #289: MSE= 14879.093223332717
Iteration #290: MSE= 14876.31146505587
Iteration #291: MSE= 14865.128027043767
Iteration #292: MSE= 14862.664259890578
Iteration #293: MSE= 14851.921722920672
Iteration #294: MSE= 14848.371201513535
Iteration #295: MSE= 14837.657635308185
Iteration #296: MSE= 14834.687399575452
Iteration #297: MSE= 14824.461543900603
Iteration #298: MSE= 14822.245578960115
Iteration #299: MSE= 14813.143058230637
Iteration #300: MSE= 14810.731450540767
Iteration #301: MSE= 14800.878607439694
Iteration #302: MSE= 14798.081284654872
Iteration #303: MSE= 14789.181334456998
Iteration #304: MSE= 14786.955954275983
Iteration #305: MSE= 14777.487212299087
Iteration #306: MSE= 14774.968910753127
Iteration #307: MSE= 14766.439804434205
Iteration #308: MSE= 14764.816207856351
Iteration #309: MSE= 14755.613942920982
Iteration #310: MSE= 14753.443151308926
Iteration #311: MSE= 14743.453067650824
Iteration #312: MSE= 14741.725410652674
Iteration #313: MSE= 14732.779952376051
Iteration #314: MSE= 14731.080001107563
Iteration #315: MSE= 14722.398336351967
Iteration #316: MSE= 14721.025982619012
Iteration #317: MSE= 14711.434933777673
Iteration #318: MSE= 14709.914048794442
Iteration #319: MSE= 14701.252249370285
Iteration #320: MSE= 14698.96975656556
Iteration #321: MSE= 14690.94520984386
Iteration #322: MSE= 14689.068340390802
Iteration #323: MSE= 14680.436206181841
Iteration #324: MSE= 14678.480401190543
Iteration #325: MSE= 14670.410456267866
Iteration #326: MSE= 14668.717030462152
Iteration #327: MSE= 14660.75437298348
Iteration #328: MSE= 14658.663327661468
Iteration #329: MSE= 14650.071236345375
Iteration #330: MSE= 14648.503131611511
Iteration #331: MSE= 14640.465960606854
Iteration #332: MSE= 14638.009267871212
Iteration #333: MSE= 14628.857324973023
Iteration #334: MSE= 14626.67834122902
Iteration #335: MSE= 14618.747572055225
Iteration #336: MSE= 14616.973223997942
Iteration #337: MSE= 14609.171701733261
Iteration #338: MSE= 14607.713229412178
Iteration #339: MSE= 14599.842022670447
Iteration #340: MSE= 14598.417399030297
Iteration #341: MSE= 14589.940722835196
Iteration #342: MSE= 14588.354728190427
Iteration #343: MSE= 14580.244802421248
Iteration #344: MSE= 14578.968077489324
Iteration #345: MSE= 14571.210383359443
Iteration #346: MSE= 14569.686697904353
Iteration #347: MSE= 14562.869630398956
Iteration #348: MSE= 14562.34226875458
Iteration #349: MSE= 14554.180949889407
Iteration #350: MSE= 14552.60367236601
Iteration #351: MSE= 14546.96461277427
Iteration #352: MSE= 14545.52188559763
Iteration #353: MSE= 14538.46839725189
Iteration #354: MSE= 14536.901757868212
Iteration #355: MSE= 14530.263029059091
Iteration #356: MSE= 14528.484853224292
Iteration #357: MSE= 14522.455942429573
Iteration #358: MSE= 14521.884549909768
Iteration #359: MSE= 14515.535776517396
Iteration #360: MSE= 14515.110074549486
Iteration #361: MSE= 14508.979616087048
Iteration #362: MSE= 14508.151475197024
Iteration #363: MSE= 14502.530112241322
Iteration #364: MSE= 14500.367201557287
Iteration #365: MSE= 14493.339521368667
Iteration #366: MSE= 14492.531973022938
Iteration #367: MSE= 14486.29288766241
Iteration #368: MSE= 14485.472778962987
Iteration #369: MSE= 14479.261746033815
Iteration #370: MSE= 14478.479622482693
Iteration #371: MSE= 14472.502908713977
Iteration #372: MSE= 14470.998444531786
Iteration #373: MSE= 14465.047490298395
Iteration #374: MSE= 14465.387767863645 so we use the result of the previous iteration.

#collapse_hide
# We can keep going the iteration until lowest MSE

#change alpha if you like
alpha=2

i=0
while i < 10:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  print(msecalc(Ipd[:len(conf_df_pd.loc['France'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['France'])]))
  i=i+1

14459.852568613716
14459.351561067231
14453.957064700979
14453.452195254797
14447.490444532958
14446.265432120781
14439.930296478404
14439.300086843774
14433.80505894358
14433.59714741496

Code

fig = go.Figure(data=[    
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['France'].index, y=Enext),
    go.Scatter(name='Ipd=conv(deconv(Ipd))', x=inf_df.loc['France'].index, y=signal.fftconvolve(h_L, Enext, mode='full')),
    go.Bar(name='Ipd', x=inf_df.loc['France'].index, y=conf_df_pd.loc['France']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=inf_df.loc['France'].index, y=lowpass(conf_df_pd.loc['France'], 0.02, 1)[0])
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{France: Actual } I_{pd} \text{ vs. convolution of deconvolution of } I_{pd}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

3. \(\beta\) and \(R\) from \(E_{pd}\) and \(I\)

As described above:

\[R = \frac{E_{pd}}{\gamma~I}\]

Code

# Calculate R
gam = 1/20.11 # As we say gamma is 1/20.11
R = Enext[:len(inf_df.loc['France'])]*(1/gam)/inf_df.loc['France']

fig = go.Figure(data=[    
    go.Scatter(name='R', x=inf_df.loc['France'].index, y=R),
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['France'].index, y=Enext),
    go.Scatter(name='Inf', x=inf_df.loc['France'].index, y=inf_df.loc['France']),
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'R',
    title={
        'text':r'$\text{France: R }$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

After March 15th we see a rapi decline in \(R\) until \(R<1\) since April 12th, however the last 2 weeks of May have seen an increase and it is now close to 1 in early June.

Germany

Use fc = 0.05

1. Checking \(h_I\)

Code

fig = go.Figure(data=[    
    go.Bar(name='Ipd', x=conf_df_pd.loc['Germany'].index, y=conf_df_pd.loc['Germany']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=conf_df_pd.loc['Germany'].index, y=lowpass(conf_df_pd.loc['Germany'], 0.05, 1)[0]),
    go.Bar(name='Rpd', x=rec_df_pd.loc['Germany'].index, y=rec_df_pd.loc['Germany']),
    go.Scatter(name='Rpd=lowpass(Rpd)', x=rec_df_pd.loc['Germany'].index, y=lowpass(rec_df_pd.loc['Germany'], 0.05, 1)[0]),
    go.Scatter(name='Rpd=conv(Ipd)', x=conf_df_pd.loc['Germany'].index, y=signal.fftconvolve(h_I, conf_df_pd.loc['Germany'], mode='full'))
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Germany: Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see the actual \(R_{pd}\) lags \(h_I[j]\circledast I_{pd}[j]\) by about 10 days.

Although they are pretty close still.

2. Estimating \(E_{pd}\) by deconvolution of \(I_{pd}\)

Code

#Settting up for deconvolution of Ipd

#regularization parameter
alpha=2

# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(conf_df_pd.loc['Germany'], 0.05, 1)[0])
Ipd[Ipd<0]=0


# Pad with last value
i=0
while i < 100:
  Ipd=np.append(Ipd, Ipd[-1])
  i=i+1

# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()

# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess

# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 10000000)
mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Germany'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Germany'])]))

itercount=0
while mse[-1] < mse[-2]:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Germany'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Germany'])]))
  print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")

Iteration #1: MSE= 4480813.327614089
Iteration #2: MSE= 1779548.6584457145
Iteration #3: MSE= 1047028.0176745327
Iteration #4: MSE= 763775.8058749898
Iteration #5: MSE= 616594.4567051543
Iteration #6: MSE= 529747.000092582
Iteration #7: MSE= 476190.4914439757
Iteration #8: MSE= 439455.98413487274
Iteration #9: MSE= 413238.4188322324
Iteration #10: MSE= 393123.1685706109
Iteration #11: MSE= 377495.96667679405
Iteration #12: MSE= 364719.9359760672
Iteration #13: MSE= 354393.91268630454
Iteration #14: MSE= 345599.6601037873
Iteration #15: MSE= 338346.69244214304
Iteration #16: MSE= 331946.96865472384
Iteration #17: MSE= 326628.71291526343
Iteration #18: MSE= 321807.5307918656
Iteration #19: MSE= 317817.5391580721
Iteration #20: MSE= 314080.2491777902
Iteration #21: MSE= 311017.61641710735
Iteration #22: MSE= 308056.9970819047
Iteration #23: MSE= 305677.4624175552
Iteration #24: MSE= 303277.81373376475
Iteration #25: MSE= 301399.156928807
Iteration #26: MSE= 299424.407404202
Iteration #27: MSE= 297923.86043155007
Iteration #28: MSE= 296276.72521180415
Iteration #29: MSE= 295073.7223713621
Iteration #30: MSE= 293668.5313263045
Iteration #31: MSE= 292701.9214976401
Iteration #32: MSE= 291486.3663743539
Iteration #33: MSE= 290701.1734301538
Iteration #34: MSE= 289654.7917024821
Iteration #35: MSE= 288995.9194548779
Iteration #36: MSE= 288056.58729253773
Iteration #37: MSE= 287539.23734438093
Iteration #38: MSE= 286709.20267507294
Iteration #39: MSE= 286278.69113009755
Iteration #40: MSE= 285518.0192174311
Iteration #41: MSE= 285165.7646299935
Iteration #42: MSE= 284483.7725251523
Iteration #43: MSE= 284191.1919826371
Iteration #44: MSE= 283555.7338814083
Iteration #45: MSE= 283332.78778147144
Iteration #46: MSE= 282755.7882248326
Iteration #47: MSE= 282569.3144878955
Iteration #48: MSE= 282028.92531054956
Iteration #49: MSE= 281873.5443270758
Iteration #50: MSE= 281368.16484932543
Iteration #51: MSE= 281260.0466810068
Iteration #52: MSE= 280777.22983139596
Iteration #53: MSE= 280702.1100166799
Iteration #54: MSE= 280245.97881903534
Iteration #55: MSE= 280179.61994918366
Iteration #56: MSE= 279744.2434719691
Iteration #57: MSE= 279711.2099280925
Iteration #58: MSE= 279298.9192786068
Iteration #59: MSE= 279292.322029075
Iteration #60: MSE= 278889.48600294883
Iteration #61: MSE= 278901.1450810871 so we use the result of the previous iteration.

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Code

# We can keep going the iteration until lowest MSE

#change alpha if you like
alpha=2

i=0
while i < 10:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  print(msecalc(Ipd[:len(conf_df_pd.loc['Germany'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Germany'])]))
  i=i+1

278505.5331361171
278534.93735467247
278160.17827383784
278200.11201937246
277842.2676195372
277879.94690998486
277532.45464883576
277593.0411346929
277256.49456925044
277320.03160332504

Code

fig = go.Figure(data=[    
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Germany'].index, y=Enext),
    go.Scatter(name='Ipd=conv(deconv(Ipd))', x=inf_df.loc['Germany'].index, y=signal.fftconvolve(h_L, Enext, mode='full')),
    go.Bar(name='Ipd', x=inf_df.loc['Germany'].index, y=conf_df_pd.loc['Germany']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=inf_df.loc['Germany'].index, y=lowpass(conf_df_pd.loc['Germany'], 0.05, 1)[0])
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Germany: Actual } I_{pd} \text{ vs. convolution of deconvolution of } I_{pd}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see that our estimate for \(\hat{E}_{pd}\) must be close to the reality of \(E_{pd}\) as \(I_{pd}\) is almost identical to \(\hat{E}_{pd} \circledast h_L\).

Of course, this holds only as long as our estimate of \(h_L\) is close to reality.

3. \(\beta\) and \(R\) from \(E_{pd}\) and \(I\)

As described above:

\[R = \frac{E_{pd}}{\gamma~I}\]

Code

# Calculate R
gam = 1/20.11 # As we say gamma is 1/20.11
R = Enext[:len(inf_df.loc['Germany'])]*(1/gam)/inf_df.loc['Germany']

fig = go.Figure(data=[    
    go.Scatter(name='R', x=inf_df.loc['Germany'].index, y=R),
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Germany'].index, y=Enext),
    go.Scatter(name='Inf', x=inf_df.loc['Germany'].index, y=inf_df.loc['Germany']),
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'R',
    title={
        'text':r'$\text{Germany: R }$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

We can see \(R<1\) since April 3rd, however the peak \(E_{pd}\) is slightly ahead of that March 23rd.

Switzerland

1. Checking \(h_I\)

Code

fig = go.Figure(data=[    
    go.Bar(name='Ipd', x=conf_df_pd.loc['Switzerland'].index, y=conf_df_pd.loc['Switzerland']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=conf_df_pd.loc['Switzerland'].index, y=lowpass(conf_df_pd.loc['Switzerland'], 0.05, 1)[0]),
    go.Bar(name='Rpd', x=rec_df_pd.loc['Switzerland'].index, y=rec_df_pd.loc['Switzerland']),
    go.Scatter(name='Rpd=lowpass(Rpd)', x=rec_df_pd.loc['Switzerland'].index, y=lowpass(rec_df_pd.loc['Switzerland'], 0.05, 1)[0]),
    go.Scatter(name='Rpd=conv(Ipd)', x=conf_df_pd.loc['Switzerland'].index, y=signal.fftconvolve(h_I, conf_df_pd.loc['Switzerland'], mode='full'))
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Switzerland: Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see the actual \(R_{pd}\) leads \(h_I[j]\circledast I_{pd}[j]\) by about 7 days.

2. Estimating \(E_{pd}\) by deconvolution of \(I_{pd}\)

Code

#Settting up for deconvolution of Ipd

#regularization parameter
alpha=2

# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(conf_df_pd.loc['Switzerland'], 0.05, 1)[0])
Ipd[Ipd<0]=0


# Pad with last value
i=0
while i < 100:
  Ipd=np.append(Ipd, Ipd[-1])
  i=i+1

# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()

# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess

# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 10000000)
mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Switzerland'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Switzerland'])]))

itercount=0
while mse[-1] < mse[-2]:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['Russia'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Switzerland'])]))
  print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")

Iteration #1: MSE= 195938.77912974023
Iteration #2: MSE= 92598.5502258455
Iteration #3: MSE= 59128.50699903844
Iteration #4: MSE= 42822.838779538215
Iteration #5: MSE= 34413.74497872791
Iteration #6: MSE= 29397.444434795136
Iteration #7: MSE= 26266.214258011416
Iteration #8: MSE= 24082.409100474433
Iteration #9: MSE= 22492.896300229964
Iteration #10: MSE= 21225.357922625102
Iteration #11: MSE= 20202.622330079314
Iteration #12: MSE= 19310.44685085919
Iteration #13: MSE= 18545.582296652286
Iteration #14: MSE= 17842.045548122835
Iteration #15: MSE= 17220.157891847626
Iteration #16: MSE= 16635.43422388673
Iteration #17: MSE= 16109.62920614112
Iteration #18: MSE= 15607.948586339993
Iteration #19: MSE= 15157.358178709917
Iteration #20: MSE= 14719.405929680112
Iteration #21: MSE= 14321.39639162879
Iteration #22: MSE= 13937.747258224457
Iteration #23: MSE= 13580.699248153194
Iteration #24: MSE= 13238.287487807323
Iteration #25: MSE= 12923.05672740925
Iteration #26: MSE= 12615.873340651639
Iteration #27: MSE= 12333.948423913023
Iteration #28: MSE= 12060.586862965463
Iteration #29: MSE= 11806.337909513775
Iteration #30: MSE= 11556.829332708503
Iteration #31: MSE= 11324.758908874524
Iteration #32: MSE= 11097.827572308348
Iteration #33: MSE= 10889.767639511578
Iteration #34: MSE= 10684.970618385503
Iteration #35: MSE= 10495.827476306413
Iteration #36: MSE= 10309.372490124444
Iteration #37: MSE= 10136.62614083432
Iteration #38: MSE= 9962.25527682353
Iteration #39: MSE= 9806.712548334388
Iteration #40: MSE= 9645.55837753218
Iteration #41: MSE= 9504.786823092121
Iteration #42: MSE= 9360.478542703051
Iteration #43: MSE= 9227.271539507645
Iteration #44: MSE= 9091.223107875738
Iteration #45: MSE= 8969.29115341354
Iteration #46: MSE= 8844.02110950773
Iteration #47: MSE= 8731.449306161778
Iteration #48: MSE= 8614.192146629004
Iteration #49: MSE= 8512.411724441054
Iteration #50: MSE= 8402.36013411898
Iteration #51: MSE= 8306.971958009977
Iteration #52: MSE= 8206.804664945614
Iteration #53: MSE= 8116.658953537234
Iteration #54: MSE= 8021.574930314654
Iteration #55: MSE= 7938.3376066661285
Iteration #56: MSE= 7847.576233211533
Iteration #57: MSE= 7770.313752856979
Iteration #58: MSE= 7684.890990888229
Iteration #59: MSE= 7614.203469219121
Iteration #60: MSE= 7531.785696940037
Iteration #61: MSE= 7465.689114837634
Iteration #62: MSE= 7387.706935147723
Iteration #63: MSE= 7327.12672350847
Iteration #64: MSE= 7253.531286391747
Iteration #65: MSE= 7197.811648480507
Iteration #66: MSE= 7128.7513445972445
Iteration #67: MSE= 7077.0580090597905
Iteration #68: MSE= 7008.087621851146
Iteration #69: MSE= 6959.073688894325
Iteration #70: MSE= 6893.160330404095
Iteration #71: MSE= 6848.244073372437
Iteration #72: MSE= 6783.542205063829
Iteration #73: MSE= 6741.731420071753
Iteration #74: MSE= 6675.396365489819
Iteration #75: MSE= 6636.0554391552005
Iteration #76: MSE= 6572.390374418994
Iteration #77: MSE= 6536.405026680437
Iteration #78: MSE= 6470.231534183589
Iteration #79: MSE= 6437.400466541315
Iteration #80: MSE= 6369.952847739624
Iteration #81: MSE= 6337.158752046109
Iteration #82: MSE= 6267.110753591126
Iteration #83: MSE= 6234.511459392571
Iteration #84: MSE= 6164.3484413572905
Iteration #85: MSE= 6135.093346091662
Iteration #86: MSE= 6064.28180181078
Iteration #87: MSE= 6035.566267853741
Iteration #88: MSE= 5962.446284991836
Iteration #89: MSE= 5934.705069962215
Iteration #90: MSE= 5860.071406975256
Iteration #91: MSE= 5833.161439680613
Iteration #92: MSE= 5757.577703215929
Iteration #93: MSE= 5734.022807025684
Iteration #94: MSE= 5656.5955806182155
Iteration #95: MSE= 5632.317638378359
Iteration #96: MSE= 5555.332007992863
Iteration #97: MSE= 5533.556372087182
Iteration #98: MSE= 5458.3969114889
Iteration #99: MSE= 5436.268504588962
Iteration #100: MSE= 5359.929024443102
Iteration #101: MSE= 5341.18226004137
Iteration #102: MSE= 5266.790948124076
Iteration #103: MSE= 5250.868203079338
Iteration #104: MSE= 5177.763832631689
Iteration #105: MSE= 5162.957706817038
Iteration #106: MSE= 5092.904915671228
Iteration #107: MSE= 5080.8933161121495
Iteration #108: MSE= 5014.749809938904
Iteration #109: MSE= 5005.743957352035
Iteration #110: MSE= 4943.5104626660495
Iteration #111: MSE= 4937.749021715691
Iteration #112: MSE= 4878.624288317557
Iteration #113: MSE= 4875.938516831257
Iteration #114: MSE= 4820.54753118064
Iteration #115: MSE= 4818.20354907714
Iteration #116: MSE= 4767.369486952111
Iteration #117: MSE= 4769.003831843811 so we use the result of the previous iteration.

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Code

# We can keep going the iteration until lowest MSE

#change alpha if you like
alpha=2

i=0
while i < 10:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  print(msecalc(Ipd[:len(conf_df_pd.loc['Switzerland'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['Switzerland'])]))
  i=i+1

4723.25082742416
4725.699669287725
4683.134053574572
4687.788070251505
4650.55132533759
4656.470026778366
4620.760015152314
4626.699797431839
4596.550489682245
4604.349708285622

Code

fig = go.Figure(data=[    
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Switzerland'].index, y=Enext),
    go.Scatter(name='Ipd=conv(deconv(Ipd))', x=inf_df.loc['Switzerland'].index, y=signal.fftconvolve(h_L, Enext, mode='full')),
    go.Bar(name='Ipd', x=inf_df.loc['Switzerland'].index, y=conf_df_pd.loc['Switzerland']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=inf_df.loc['Switzerland'].index, y=lowpass(conf_df_pd.loc['Switzerland'], 0.05, 1)[0])
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{Switzerland: Actual } I_{pd} \text{ vs. convolution of deconvolution of } I_{pd}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

We can see that our estimate for \(\hat{E}_{pd}\) must be close to the reality of \(E_{pd}\) as \(I_{pd}\) is almost identical to \(\hat{E}_{pd} \circledast h_L\).

Of course, this holds only as long as our estimate of \(h_L\) is close to reality.

3. \(\beta\) and \(R\) from \(E_{pd}\) and \(I\)

As described above:

\[R = \frac{E_{pd}}{\gamma~I}\]

#collapse_hide

# Calculate R
gam = 1/20.11 # As we say gamma is 1/20.11
R = Enext[:len(inf_df.loc['Switzerland'])]*(1/gam)/inf_df.loc['Switzerland']

fig = go.Figure(data=[    
    go.Scatter(name='R', x=inf_df.loc['Switzerland'].index, y=R),
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['Switzerland'].index, y=Enext),
    go.Scatter(name='Inf', x=inf_df.loc['Switzerland'].index, y=inf_df.loc['Switzerland']),
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'R',
    title={
        'text':r'$\text{Switzerland: R }$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

United States

1. Checking \(h_I\)

Code

fig = go.Figure(data=[    
    go.Bar(name='Ipd', x=conf_df_pd.loc['US'].index, y=conf_df_pd.loc['US']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=conf_df_pd.loc['US'].index, y=lowpass(conf_df_pd.loc['US'], 0.05, 1)[0]),
    go.Bar(name='Rpd', x=rec_df_pd.loc['US'].index, y=rec_df_pd.loc['US']),
    go.Scatter(name='Rpd=lowpass(Rpd)', x=rec_df_pd.loc['US'].index, y=lowpass(rec_df_pd.loc['US'], 0.05, 1)[0]),
    go.Scatter(name='Rpd=conv(Ipd)', x=conf_df_pd.loc['US'].index, y=signal.fftconvolve(h_I, conf_df_pd.loc['US'], mode='full'))
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{US: Actual } R_{pd} \text{ vs. } h_I[j]\circledast I_{pd}[j]$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Data is not very good.

2. Estimating \(E_{pd}\) by deconvolution of \(I_{pd}\)

Code

#Settting up for deconvolution of Ipd

#regularization parameter
alpha=2

# Setup up the resultant Ipd we want to compare our guess with
Ipd=np.floor(lowpass(conf_df_pd.loc['US'], 0.05, 1)[0])
Ipd[Ipd<0]=0


# Pad with last value
i=0
while i < 100:
  Ipd=np.append(Ipd, Ipd[-1])
  i=i+1

# Find delay caused by h_L
delay=Ipd.argmax()-signal.fftconvolve(Ipd, h_L, mode='full').argmax()

if (np.abs(delay)>20):
  delay = -15

# We want initial guess to simply be the result of the convolution delayed
initial_guess=np.roll(Ipd,delay)
Enext = initial_guess

# AN array to record MSE between result we want and our iterated guess
mse=np.array([])
mse=np.append(mse, 10000000)
mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['US'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['US'])]))

itercount=0
while mse[-1] < mse[-2]:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  mse=np.append(mse, msecalc(Ipd[:len(conf_df_pd.loc['US'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['US'])]))
  print("Iteration #" + str(itercount) +": MSE= "+str(mse[itercount]))
print("Iteration #" + str(itercount+1) +": MSE= "+str(mse[-1])+" so we use the result of the previous iteration.")

Iteration #1: MSE= 47992077.600439996 so we use the result of the previous iteration.

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

Code

# We can keep going the iteration until lowest MSE

#change alpha if you like
alpha=2

i=0
while i < 10:
  itercount=itercount+1
  Enext=iter_deconv(alpha, h_L, Enext, delay, Ipd)
  print(msecalc(Ipd[:len(conf_df_pd.loc['US'])], signal.fftconvolve(h_L, Enext, mode='full')[:len(conf_df_pd.loc['US'])]))
  i=i+1

22255749.4520845
13474646.870988503
9247852.125484189
6845961.888525094
5353233.801469828
4373948.068754885
3688854.4786455603
3193666.3164928
2820280.9677202436
2533321.732436159

Code

fig = go.Figure(data=[    
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['US'].index, y=Enext),
    go.Scatter(name='Ipd=conv(deconv(Ipd))', x=inf_df.loc['US'].index, y=signal.fftconvolve(h_L, Enext, mode='full')),
    go.Bar(name='Ipd', x=inf_df.loc['US'].index, y=conf_df_pd.loc['US']),
    go.Scatter(name='Ipd=lowpass(Ipd)', x=inf_df.loc['US'].index, y=lowpass(conf_df_pd.loc['US'], 0.05, 1)[0])
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'Count',
    title={
        'text':r'$\text{US: Actual } I_{pd} \text{ vs. convolution of deconvolution of } I_{pd}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

/tmp/ipykernel_17425/3295257807.py:10: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

3. \(\beta\) and \(R\) from \(E_{pd}\) and \(I\)

As described above:

\[R = \frac{E_{pd}}{\gamma~I}\]

Code

# Calculate R
gam = 1/20.11 # As we say gamma is 1/20.11
R = Enext[:len(inf_df.loc['US'])]*(1/gam)/inf_df.loc['US']

fig = go.Figure(data=[    
    go.Scatter(name='R', x=inf_df.loc['US'].index, y=R),
    go.Scatter(name='Epd=deconv(Ipd)', x=inf_df.loc['US'].index, y=Enext),
    go.Scatter(name='Inf', x=inf_df.loc['US'].index, y=inf_df.loc['US']),
])

fig.update_layout(
    
    xaxis_title = 'Day',
    yaxis_title = 'R',
    title={
        'text':r'$\text{US: R }$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

After March 15th we see a rapi decline in \(R\) until \(R<1\) since April 12th, however the last 2 weeks of May have seen an increase and it is now close to 1 in early June.

We can see \(R\) declined rapidly after March 12th to \(R<1\) on May 30th, but has since grown back to close to 1.