Source: https://github.com/llorracc/EpiExp/blob/master/SIR_Ndlib.ipynb
This notebook is based on the model in Shiller and Pound (1989)
Author: Tao Wang and Chris Carroll
Date: September, 2021
This notebook uses NDlib to simulate the variant described in Carroll and Wang (2021) of the epidemiological model of expectations of Shiller and Pound (1989).
We exploit the fact that for a large population, a random graph
of the kind that the underlying NetworkX library handles is a good approximation to the discrete-time SIR model, when appropriately parameterized.
$N$: A large population of investors, divided into three "compartments"
$\beta$: what epidemiologists call the infection rate
random mixing
and the disease transmission upon encounters between $I$ and $S$\begin{equation} \beta =\underbrace{\chi}_{\text{contact rate}} \times \underbrace{\tau}_{\text{transmission probability}} \end{equation}
$\gamma$: recovery rate
Putting these together, the change in the population of infected persons is given by
\begin{equation} I_{t+1}- I_t = \beta \frac{S_t}{N}I_t - \gamma I_t \end{equation}The first term on the right hand side above captures the number of people who "flow" from 'compartment $S$' to 'compartment $I$', which is proportional to the infection rate $\beta$, the fraction of people who are susceptible $\frac{S_t}{N}$, and the number of the infected $I_t$. The second term captures the number of people who "flow" from $I$ to $R$ (See the flow diagram).
Our treatment makes two modifications to the original model as described in Shiller and Pound (1989). First, for consistency with the software toolkit used below, we rewrite their originally continuous-time model in a discrete time form. Second, the original paper described an additional stochastic shock to the change in $I_t$ meant to capture a potential "change in the 'source' of the infection or the nature of the contagion." Because that shock was not actually used for any results in the paper, we neglect it.
An analytical formula for the limiting size of the compartment $R$ exists under following assumptions. (See the Wikipedia Page)
See the Appendix for the detailed derivations.
import ndlib # 'network dynamics library'
import networkx as nx # built on top of networkx toolkit
import matplotlib.pyplot as plt
import ndlib.models.epidemics as ep
import ndlib.models.ModelConfig as mc
from ndlib.viz.mpl.DiffusionTrend import DiffusionTrend
from ndlib.viz.mpl.DiffusionPrevalence import DiffusionPrevalence
from ndlib.viz.bokeh.MultiPlot import MultiPlot
%matplotlib inline
from SIR import SIR, solve_r_ss
import numpy as np
no display found. Using non-interactive Agg backend
## plot the model
model = nx.DiGraph()
model.add_edges_from([("S", "I"), ("I", "R")])
plt.figure(figsize=(12, 3))
pos = {'S':(0, 0.0), 'I':(10, 0.0), 'R':(20, 0.0)}
nx.draw_networkx(model,
pos = pos,
arrows=True,
node_size = 4000,
node_shape ='s',
node_color= 'gainsboro',
alpha = 0.5,
font_weight = 40,
font_size = 35,
arrowsize = 30,
arrowstyle = 'fancy',
)
nx.draw_networkx_edge_labels(model,
pos,
edge_labels={('S','I'):r'$\beta \frac{S_t}{N}I_t$ ',\
('I','R'):r'$\gamma I_t$'},
font_color='black',
font_size = 35)
# See if you can add a bit of space between the equations and the boxes
# Add time subscripts in the flow diagram
plt.savefig("./draft/chapter/figures/flow_diagram.png",format="PNG")
To illustrate the model’s implications under such a setting, we parameterize the model with four such combinations of parameter values taken from Shiller and Pound (1989), characterizing two different kinds of investors and two categories of stocks.
We convert all the continuous-time rates into discrete-time and from annual to weekly frequency. For instance, the recovery rate estimated from the decaying pattern of the time spent on studying a given stock for INSRPI is $g=1.39$ (a half-life of $ln(2)/g=0.50$ years). In discrete-time and at weekly frequency, this is equivalent to a probability of recovery $\gamma = 1-\exp^{-g/52} =0.02$. For the infection rate, under the assumption made by the paper that the fraction of susceptible is close to $1$ despite being time-varying, the estimated median infection rate of INSRPI is $b = 2.02$. It is converted to a weekly probability of $\beta = 1-\exp^{-b/52}=0.038$. We set the initial fraction of the infected to be $1$ percent.
## a function that converts continuous rate to discrete-time weekly rates
def ct2dtw(rate):
return 1-np.exp(-rate/52)
## Shiller and Pound estimates
types = ['INSRAND','INSRPI','INDRAND','INDRPI']
gamma_continuous_ests = [1.84,1.39,3.72,1.73] ## directly take from the paper (page 58)
beta_continuous_ests = [20.66,7.35,12.74,3.71] ## directly take from the paper (page 58)
beta_sir_list = [round(ct2dtw(beta),2) for beta in beta_continuous_ests]
gamma_list = [round(ct2dtw(gamma),2) for gamma in gamma_continuous_ests]
The toolkit we are using, NDlib, implements a SIR model residing on a given simulated network from another scientific library NetworkX.
Note that the network-based SIR model does not take the infection rate $\beta$, as a catch-all parameter. Instead, it let the users configure the transmission probability $\tau$ (what the program calls $\beta$, which is not the $\beta$ defined above), and the network structure that affects the contact rates $\chi$, separately. In an N-sized random graph a la Erdos and Renyi (1960) with a connection probability $p$, we have the following relationship
\begin{equation} \chi = \underbrace{(N-1)p}_{\text{average degree}}. \end{equation}Hence, we can use the following equation to compute the underlying transmission probability $\tau$ corresponding to the $\beta$ in a random graph parameterized by $N$ and $p$.
\begin{equation} \tau = \frac{\beta}{\chi} = \frac{\beta}{(N-1)p} \end{equation}Therefore, we write the following function to compute the value of $\tau$ corresponding to a value of $\beta$ we intend to use, depending on the contact rate $\chi$.
## beta in NDlib is not really beta in standard SIR
## the infection probability tau = beta in SIR / χ, where χ is the number of contacts
## χ = = average degree of a random graph = (N-1) x p
def beta2tau(beta_sir,
χ):
return beta_sir/χ
## create the population residing at a random graph
random_prb = 0.1
population = 8000
degree = random_prb*(population-1)
χ = degree ## number of contacts per period per person is the degree in a random graph
## get the corresponding τ values for the weekly infection rate β
τ_list = [beta2tau(beta_sir,χ) for beta_sir in beta_sir_list]
random_mix = nx.erdos_renyi_graph(population,random_prb) # agents are randomly connected
# SIR model to simulate the spread of economic opinions
sir = ep.SIRModel(random_mix)
sir.available_statuses
{'Susceptible': 0, 'Infected': 1, 'Removed': 2}
## simulation periods
iters = 208
## prepare the simulation
trends_list =[]
## simulation
nb = len(types)
for x in range(nb):
sir = ep.SIRModel(random_mix)
cfg_sir = mc.Configuration()
cfg_sir.add_model_parameter('beta', τ_list[x]) # transmission probability = infection rate/nb of contacts
cfg_sir.add_model_parameter('gamma', gamma_list[x]) # removal rate
cfg_sir.add_model_parameter("fraction_infected", 0.01)
sir.set_initial_status(cfg_sir)
## simulate the sir model for a particular configuration
iterations = sir.iteration_bunch(iters, node_status=True)
trends_sir = sir.build_trends(iterations)
trends_list.append(trends_sir)
100%|█████████████████████████████████████████| 208/208 [00:25<00:00, 8.22it/s] 100%|█████████████████████████████████████████| 208/208 [00:35<00:00, 5.81it/s] 100%|█████████████████████████████████████████| 208/208 [00:18<00:00, 11.29it/s] 100%|█████████████████████████████████████████| 208/208 [00:42<00:00, 4.94it/s]
## some parameters to testing
beta = 0.1
gamma = 0.03
i0 = 0.01
s0 = 1-i0
r0 = 0.0
x0 = (s0,i0,r0)
T = iters
times = range(T)
## compare the paths from Ndlib and SIR
lw = 5
lbsize = 15
## plot
fig, axs = plt.subplots(4,1,
figsize=(13, 22),
facecolor='w',
edgecolor='k')
fig.subplots_adjust(hspace = 0.3, wspace=.1)
axs = axs.ravel()
nb = len(types)
nbw1y = 52
for x in range(nb):
beta = beta_sir_list[x]
gamma = gamma_list[x]
s,i,r = SIR(beta, # infection rate
gamma, # recovery rate
T,
x0)
title = types[x]+': '+r'$\beta={}$'.format(round(beta,2))+','+r'$\gamma ={}$'.format(round(gamma,2))
axs[x].set_title(title,
fontsize=lbsize+5)
r_ss = solve_r_ss(beta,gamma)
axs[x].hlines(r_ss,
0.0,
T,
color='r',
linestyle ='solid',
lw = 1,
label=r'$R_{+\infty}$')
### generated from NDlib
s_nd = np.array(trends_list[x][0]['trends']['node_count'][0])/population
i_nd = np.array(trends_list[x][0]['trends']['node_count'][1])/population
r_nd = np.array(trends_list[x][0]['trends']['node_count'][2])/population
axs[x].plot(times,s_nd,
'--',
lw=lw,
label=r'$S_t$')
axs[x].plot(times,i_nd,
'-.',
lw=lw,
label=r'$I_t$ ')
axs[x].plot(times,r_nd,
'-',
lw=lw,
label=r'$R_t$ ')
### generated from SIR
#axs[x].plot(times,s,'-',lw=lw,label='S')
#axs[x].plot(times,i,'--',lw=lw,label='I')
#axs[x].plot(times,r,'-.',lw=lw,label='R')
axs[x].set_xlim(0.0,T)
axs[0].legend(loc=0,
prop={'size': 22})
axs[x].tick_params(axis='x', labelsize=lbsize)
axs[x].tick_params(axis='y', labelsize=lbsize)
plt.savefig("./draft/chapter/figures/sir_simulate.png",format="PNG")
An analytical formula for the limiting size of the compartment $R$ exists under following assumptions. (See the Wikipedia Page)
To obtain the result, we use lower case letter $s_t$, $i_t$, and $r_t$ to represent the fractions of each compartment.
\begin{equation} \begin{split} & s_{t+1} - s_{t} = - \beta s_t i_t \\ & i_{t+1}- i_{t} = \beta s_t i_t - \gamma i_t \\ & r_{t+1} - r_t = \gamma i_t \end{split} \end{equation}Divide the first equation above by $s_t$ into both sides and further rearrange it.
\begin{equation} \begin{split} &\Delta ln(s_{t+1}) \equiv \frac{s_{t+1} - s_{t}}{s_t} = - \beta i_t \\ & \rightarrow i_t = -\frac{\Delta ln(s_{t+1})}{\beta} \end{split} \end{equation}Plug it in the third equation above, we get
\begin{equation} \begin{split} &r_{t+1}- r_t = - \frac{\gamma}{\beta} (ln(s_{t+1})-ln(s_{t})) \\ \end{split} \end{equation}The equation above holds for any choice of $t$ and $t+\Delta t$, hence we have
\begin{equation} \begin{split} &r_{t+\Delta t}- r_t = - \frac{\gamma}{\beta} (ln(s_{t+\Delta t})-ln(s_{t})) \\ \end{split} \end{equation}Assuming the initial size of $i_0$ to be infinitely small, hence, $s_0=1$. Plug these $t=0$ values $s_0=1$ and $r_0 =0$, we obtain
\begin{equation} \begin{split} &r_{\Delta t} = \frac{\gamma}{\beta} ln(s_{\Delta t}) \\ \end{split} \end{equation}Taking the limit of the $\Delta t$ to $+ \infty$, and using the fact that in long run the fraction of the infection $i_{+\infty}=0$, and $1+s_{+\infty}=r_{+\infty}$, we arrive at the folowing equation.
\begin{equation} \begin{split} &r_{+\infty} = - \frac{\gamma}{\beta} ln(1-r_{+\infty}) \\ &\leftrightarrow e^{-\frac{\beta}{\gamma} r_{+\infty}} = 1-r_{+\infty} \end{split} \end{equation}$r_{+\infty}$ and $s_{+\infty}$ can be solved.
Notice that the limiting size does not depend on the specific values of $\beta$ and $\gamma$ but its ratio.