Following my last post, the SIR model has been completely and totally wrong. Answers follow from assumptions. It assumes a constant reproduction rate, and the virus peters out when sick people run in to recovered and immune people. That's not what's happening -- people responded by lowering the contact rate, long before we ran in to herd immunity.
I speculated last time about a model in which people respond to the severity of the disease by reducing contacts. Let's do it. (Warning: this post uses MathJax to show equations. It may not work on all devices.)
I modify the SIR model as presented by Chad Jones and JesúsFernández-Villaverde: \begin{align*} \Delta S_{t+1} & =-\beta S_{t}I_{t}/N\\ \Delta I_{t+1} & =\beta S_{t}I_{t}/N-\gamma I_{t}\\ \Delta R_{t+1} & =\gamma I_{t}-\theta R_{t}\\ \Delta D_{t+1} & =\delta\theta R_{t}\\ \Delta C_{t+1} & =(1-\delta)\theta R_{t}% \end{align*} S = susceptible, I = infected (and infectious), R resolving, i.e. sick but not infectious, D = dead, C = recovered and immune, N = population. The lags give the model momentum. Lowering the reproduction rate does not immediately stop the disease. The model uses exponential decays rather than fixed lags to capture timing. \(\beta\) is the number of contacts per day. A susceptible person meets \(\beta\) people per day. \(I/N\) of them are infected, so \(\beta S_{t}I_{t}/N\) become infected each day. We parameterize \(\beta\) in terms of the reproduction rate \(R_{0}\), \[ R_{0}=\beta/\gamma \] The number of infections from one sick person = number of contacts per day times the number of days contacts are infectious (on average).
The standard SIR model uses a constant \(\beta\) and hence a constant \(R_{0}\). The disease grows exponentially, then becomes limited by the declining number of susceptible people in the population. Each infected person runs in to recovered people, not susceptible people. The whole point is, that did not happen. We lowered \(\beta\) instead.
I model the evolution of \(\beta\) behaviorally. First, suppose people reduce their contacts in proportion to the chance of getting the disease. As people see more infectious people around, the danger of getting infected rises. They reduce their contacts proportionally to the number of infectious people. \[ \log(\beta_{t})=\log\beta_{0}-\alpha_I I_{t}/N_{t}. \] This function could also model a policy response.
However, due to the lack of testing we don't really know how many people are infectious at any time. So as a second model, suppose instead people or policy reduce contacts according to the current death rate, \[ \log(\beta_{t})=\log\beta_{0}-\alpha_D \Delta D_{t}/N. \] \(\beta\) is a rate, how many people do you bump in to per day. I use the log because it can't be negative. The log also captures the idea that early declines in \(\beta\) are easy, by eliminating superspreading activities. Later declines in \(\beta\) are more costly.
I use Chad and Jesús numbers, \(\gamma=0.2\) or 5 days of infectiousness on average, \(\theta=0.1\) implying 10 more days on average with the disease before it resolves, \(\delta=0.08\) \((0.8\%)\) death rate. They parameterize and estimate \(\beta\) \(\ \)via \(R_{0}=\beta/\gamma\). I take the original \(R_{0}% =5\), which is typical of their estimates, and implies \(\beta_{0}=\gamma R_{0}=1\). They estimate \(R_{0}^{\ast}=0.5\) so \(\beta^{\ast}=0.1\), which I will use to calibrate \(\alpha\). New York peaked at 90 deaths per million, but we will see the dynamics overshoot. So I'll pick \(\alpha\) in that case so that \(\beta=\beta^{\ast}\) at 50 daily deaths per milllion triggers \(R_{0}% =0.5\), i.e. \(\alpha_D\) solves \[ \log\left( 0.1\right) =\log\left( 1\right) -\alpha_D\times50/10^{6}. \] The death rate is about 1%, so I calibrate the infection model so that \(R_{0}=R_{0}^{\ast}=0.5\) at an infection rate of 5000 per million or 0.5%. \(\alpha_I\) solves \[ \log(0.1)=\log\left( 1\right) -\alpha_I \times 5000 / 10^{6}. \]
Here is my assumed reproduction rate as a function of deaths per million. The red dot is the calibration point: at 50 deaths per day, people and policy will drive the reproduction rate to 0.5. The red dashed line is a much more aggressive response, which I'll investigate later.
The standard SIR model
As background, here is a simulation of the standard SIR model with these numbers, and a constant \(\beta=1\) meaning \(R_0=5\).
I start at day 1 with a single infected person. The virus grows exponentially. The number infected peaks at about half the population. Around day 25 however, herd immunity starts to kick in. The number infected peaks. Sick people (resolving) peaks a bit later. The pandemic goes away almost as quickly as it came and it's over after two months. With \(R_0=5\) everyone gets it and 0.8% or 8000 people die.
This is the nightmare scenario presented to policy makers in February and caused the economic shutdowns. It is completely wrong -- it's not what happened anywhere.
The behavioral SIR model
Here is the simulation of the behavioral SIR model, in which people (or policy) reacts by lowering the contact rate in response to the number infected.
The vertical scale is different. Only about 4000 people get infected here, not 1 million! The pandemic gets going with the same exponential speed (blue line), but now once infections get up to 1000 per million we see the sharp reduction in the reproduction rate (dashed black line).
This is a lot more like what we saw! A rapid rise, to a plateau, with a much more sensible set of numbers. That's the good news. The bad news is that it goes on and on and on. The minute infections decline people slack off just enough to get it going again. Responding to infections, even though there is a lag, produces very stable dynamics.
The reproduction rate asymptotes to \(R_0=1\). This is both the good news and the bad news outlined in my last post. It doesn't get worse with second waves. But it doesn't get better either.
That result not at all related to the calibration. The reproduction rate always asymptotes to one in this model, with a steady number of infections and a steady number of deaths per day, until finally after years and years we get herd immunity and all efforts to reduce contacts are turned off.
Here is the same simulation with the much stronger response, \( alpha\) is raised by a factor of 5 to the red dashed line in my first graph. No, it's not the same graph. Notice the vertical scale. This response is much less tolerant of infections, so the overall rate of infection is much lower. But the path is exactly the same.
Being more forceful does not change the reproduction rate, which still asymptotes to one. We just trundle along with much lower infections and daily death rates.
This comparison makes nice sense of what we see, per the last post. Far different regimes give rise to essentially the same dynamics, but some at much higher and some at much lower levels.
Technology offers some hope. What happens if the costs of reducing \(\beta\) become lower over time, so people can slowly become more careful while also letting the economy grow? Widespread individual testing and tracing, for example, are ways of distancing that are less costly. In this case, we steadily move from the second to last graph to the last graph. To model that, I let \( \alpha \) vary over time, growing by a factor of 2 from time 0 to time 100,
This accounts for a plateau with a slow tail. The actual reproduction rate still is close to one, but it's just enough below one to gently let the virus decay.
Deaths and information
A big objection: here I keyed behavior to the infection rate -- people are more careful the more infected people are around. But we don't see how many people are infected. We do see deaths.
Here is the simulation when people respond to the death rate rather than the infection rate
The positive feedback does eventually control the pandemic. Each wave is smaller. And this model also trends to \( R_0 = 1\) by the same mechanism. It just takes a lot of wiggles to get there.
Information, rational exceptions, and externalitities
The contrast between the first and second graphs gives a quick policy suggestion: good information on how many people are infected in one's local area would be really helpful to avoid waves of infections.If we had just enough random testing to know how many people are infected in our local area, people and officials could follow the top graphs not the bottom graph. It is not expensive. In the model, one can back out the number of people infected from the increase in the number "resolving." The rate of hospital admissions might be a widely publicized number now available that could be a very good guess.
Widespread, available (no protocol, no prescription, just go get it, free market) testing would radically reduce the economic costs of social distancing, and end this fast. (The optimal \(\alpha\) would rise by orders of magnitude. Yes, you point to externality, why should I test myself. But you ignore economic and social demands. If such testing is available, it's really easy for customers to demand you show your test. Paul Romer is right.
Of course now we get to the delicate question of public vs. private incentives. My first model seems like a reasonable guess of how people will behave -- take actions to be careful the greater my chance of getting sick is by going out.
We want, naturally, a dynamic model in which people's actions incorporate an understanding of the dynamics. In that vein, the latter graph seems unduly pessimistic. People are pretty smart and they know that the death rate is high when the danger of going out has passed. Thus, one may well expect them to foresee the dynamics, be careful when the death rate is increasing, and slack off when it is decreasing. More generally, in this deterministic model, you can back out what the state of all the variables is if you observe one of them. Thus, the rational expectations equilibrium of this model if people want to react to the number of infections is the first one, even if they can't see infections. They can back infections out of the death data. That may be too much to hope for, but reality is likely in between.
Being careful has an externality, of course, so people following a private optimum of costly but careful behavior vs. getting sick is not necessarily the social optimum. Most economists jump quickly from this observation to calibrated time-varying lockdown policies to try to control \(\beta\). But let us not forget the other side of that coin: The public policy tools are sledgehammers, which do a poor job of controlling interactions \(\beta\) at reasonable economic cost. Really, what we have are at best exhortations to be careful in the details of daily life, plus extremely expensive business shutdowns. As a concrete example, in the model one may be tempted to advocate that officials lie about the number of infected to get people to be more careful than they would be privately. But once a lie is found out, nobody believes anything anymore, and the next step is China. I still think timely and accurate information is better.
To do list
Some more thought on functional form would be useful. Do we have any data or other ways of measuring how people behave?
Obviously a real economic model would derive these behavioral responses from a maximization problem, and consider the tradeoff between more distancing \( \beta\) and economic costs.
Optimal policy may differ most from individual behavior in the dynamics. It is not worth it to an individual to be careful early when there are few sick people around, but policy considers the effect of you getting sick on everyone who gets it from you. Optimal control of \(\beta\) beckons. But to be realistic we must include the fact that public control of \(\beta\) against private wishes will be much less efficient.
On the other papers: Chad and Jesús model social distancing, whether voluntary or by policy, via a deterministic and permanent exponential decay from a state of nature \(\beta_{0}\) to a new lower value \(\beta^{\ast}\) over a period \[ \beta_{t}=\beta_{0}e^{-\lambda t}+\beta^{\ast}(1-e^{-\lambda t}). \] The point here is to realize there is feedback, and both people and policy behavior respond to facts. Their paper fits the data so far beautifully. My goal is to think about what happens next. That \( \beta\) just sits at \(beta^\ast\) as it does in their model, seems unrealistic because people aren't going to keep distancing voluntarily or involuntarily.
Eichenbaum, Rebelo and Trabandt, have an economic model of \(\beta\). People work and shop less when they are more afraid of getting sick. But the tradeoff is not very attractive. Here is their model (solid) vs. the basic SIR model (dash)
Update:
1)Equilibrium social distancing by Flavio Toxvaerd is a simple economic model with endogenous social distancing. It also produces plateaus when people choose to be safer
2) Economists vs. epidemiologists has a long history. Economists point out that disease transmission is not a biological constant, but varies with human behavior. And human behavior varies predictably in response to incentives (and information). Many beautiful facts and stories on this point are collected in Tomas Phillipson and Richard Posner's book, Private Choices and Public Health: The AIDS Epidemic in an Economic Perspective. For example, AIDS patients in clinical trials would mix their medicines together. Half of the drug for sure is better than a 50 50 chance of nothing. Thanks to a correspondent for the reminder.
3) A Multi-Risk SIR Model with Optimally Targeted Lockdown by Daron Acemoglu, Victor Chernozhukov Michael Whinston and Ivan Werning just came out. I haven't read it yet, but it is an obvious addition to the stack for people working on epidemiology models with economic incentives. It has diverse populations and transmission mechanisms, which have long struck me as a key insight. We are not all average, and that really matters here.
4) From the comments, a many-authored paper arguing that herd immunity may be much lower. Essentially the super spreaders are more likely to get the disease, so they are more likely to be immune first. "Super spreader" includes people in nursing homes, emergency room technicians, bus drivers, etc., not just jet setting partiers.
5) I missedMacroeconomic Dynamics and Reallocation in an Epidemic by Dirk Krueger Harald Uhlig and Taojun Xie
...we distinguish goods by their degree to which they can be consumed at home rather than in a social (and thus possibly contagious) context. We demonstrate that, within the model the “Swedish solution” of letting the epidemic play out without government intervention and allowing agents to shift their sectoral behavior on their own can lead to a substantial mitigation of the economic and human costs of the COVID-19 crisis, avoiding more than 80 of the decline in output and of number of deaths within one year, compared to a model in which sectors are assumed to be homogeneous. For different parameter configurations that capture the additional social distancing and hygiene activities individuals might engage in voluntarily, we show that infections may decline entirely on their own, simply due to the individually rational re-allocation of economic activity: the curve not only just flattens, it gets reversed.6) A behavioral SIR model YouTube talk from Lones Smith
7) Systematic biases in disease forecasting – The role of behavior change by Ceyhun Eksina Keith Paarpornb Joshua S. Weitzcde
...during real-world outbreaks, individuals may modify their behavior and take preventative steps to reduce infection risk. ... we evaluate this hypothesis by comparing the dynamics arising from a simple SIR epidemic model with those from a modified SIR model in which individuals reduce contacts as a function of the current or cumulative number of cases.Thanks to Andy Atkeson for the tip. See hisA note on the economic impact of coronavirus and Talk at NBER
8)...I'm sure more updates will follow.
Code
Here is the Matlab code for my plots
close all
clear all
gam = 0.2;
thet = 0.1;
delt = 0.008;
R0 = 5;
alphaD = (0 - log(0.1))*1E6/50;
alphaI = (0 - log(0.1))*1E4/50;
beta0 = R0*gam;
N = 1E6;
T = 100;
% plot beta
figure
drate = (0:100)'/1E6;
betat = exp(log(beta0) - alphaD*drate);
betat1 = exp(log(beta0) - 5*alphaD*drate);
plot(drate*1E6, betat/gam,'linewidth',2);
hold on
plot(drate*1E6, betat1/gam,'--','linewidth',2);
hold on
plot(50, 0.5, 'o','markerfacecolor','r')
legend('Original','\alpha multiplied by 5','location','best')
xlabel('Daily deaths/million');
ylabel('Reproduction rate R0 = \beta / \gamma')
axis([0 80 0 5]);
print -dpng betafig.png
% standard SIR model
S = zeros(T,1);
I = S;
R = S;
D = S;
C = S;
S(1) = N-1;
I(1) = 1;
for t = 1:T-1;
S(t+1) = S(t) -beta0*S(t)*I(t)/N;
I(t+1) = I(t) + beta0*S(t)*I(t)/N - gam*I(t);
R(t+1) = R(t) + gam*I(t)-thet*R(t);
D(t+1) = D(t) + delt*thet*R(t);
C(t+1) = C(t) + (1-delt)*thet*R(t);
end;
figure;
plot((1:T)',[S I R 100*D C]/1E6, 'linewidth',2);
legend('Susceptible','Infected','Resolving','100 x Dead','ReCovered','location','best')
xlabel('Days')
ylabel('Millions');
axis([10 70 0 1]);
title('SIR model, constant R0 = 5');
print -dpng std_sir.png
figure;
plot((2:T)',[ I(2:T) -(S(2:T)-S(1:T-1)) 100*(D(2:T)-D(1:T-1))]/1E6, 'linewidth',2);
legend('Infected','New Infections','100 x New Dead','location','best')
xlabel('Days')
ylabel('Millions');
axis([10 70 0 0.6]);
title('SIR model, constant R0 = 5');
print -dpng std_sir_diffs.png
disp('last s i r d');
disp([S(T) I(T) R(T) D(T)]);
% my model, response to infections
S = zeros(T,1);
I = S;
R = S;
D = S;
C = S;
betat = S;
S(1) = N-1;
I(1) = 1;
I(1) = 1;
S1 = S;
I1 = I;
R1 = R;
D1 = D;
C1 = C;
betat1 = betat;
S2 = S;
I2 = I;
R2 = R;
D2 = D;
C2 = C;
betat2 = betat;
for t = 1:T-1;
betat(t) = exp(log(beta0) - alphaI*((I(t))/N));
S(t+1) = S(t) -betat(t)*S(t)*I(t)/N;
I(t+1) = I(t) + betat(t)*S(t)*I(t)/N - gam*I(t);
R(t+1) = R(t) + gam*I(t)-thet*R(t);
D(t+1) = D(t) + delt*thet*R(t);
C(t+1) = C(t) + (1-delt)*thet*R(t);
betat1(t) = exp(log(beta0) - 5*alphaI*((I1(t))/N));
S1(t+1) = S1(t) -betat1(t)*S1(t)*I1(t)/N;
I1(t+1) = I1(t) + betat1(t)*S1(t)*I1(t)/N - gam*I1(t);
R1(t+1) = R1(t) + gam*I1(t)-thet*R1(t);
D1(t+1) = D1(t) + delt*thet*R1(t);
C1(t+1) = C1(t) + (1-delt)*thet*R1(t);
betat2(t) = exp(log(beta0) - (1+1*t/T)*alphaI*((I2(t))/N));
S2(t+1) = S2(t) -betat2(t)*S2(t)*I2(t)/N;
I2(t+1) = I2(t) + betat2(t)*S2(t)*I2(t)/N - gam*I2(t);
R2(t+1) = R2(t) + gam*I2(t)-thet*R2(t);
D2(t+1) = D2(t) + delt*thet*R2(t);
C2(t+1) = C2(t) + (1-delt)*thet*R2(t);
end;
figure;
%yyaxis left
plot((2:T)',[ I(2:T) 100*(D(2:T)-D(1:T-1)) ],'linewidth',2)
xlabel('Days')
ylabel('People');
axis([0 70 0 inf]);
yyaxis right
plot((1:T)',betat/gam,'--k','linewidth',2);
ylabel('Reproduction rate R0','color','k')
axis([0 70 0 inf]);
legend('Infected','100 x Deaths/day','R0 (right scale)','location','best')
title('BSIR model, R0 varies with infection rate');
print -dpng std_sir_aI.png
figure;
%yyaxis left
plot((2:T)',[ I1(2:T) 100*(D1(2:T)-D1(1:T-1)) ],'linewidth',2)
xlabel('Days')
ylabel('People');
axis([0 70 0 inf]);
yyaxis right
plot((1:T)',betat1/gam,'--k','linewidth',2);
ylabel('Reproduction rate R0','color','k')
axis([0 70 0 inf]);
legend('Infected','100 x Deaths/day','R0 (right scale)','location','best')
title('BSIR model, R0 varies with infection rate, higher \alpha');
print -dpng std_sir_aI1.png
figure;
%yyaxis left
plot((2:T)',[ I2(2:T) 100*(D2(2:T)-D2(1:T-1)) ],'linewidth',2)
xlabel('Days')
ylabel('People');
axis([0 70 0 inf]);
yyaxis right
plot((1:T)',betat2/gam,'--k','linewidth',2);
ylabel('Reproduction rate R0','color','k')
axis([0 70 0 inf]);
legend('Infected','100 x Deaths/day','R0 (right scale)','location','best')
title('BSIR model, R0 varies with infection rate, \alpha increases over time');
print -dpng std_sir_aI2.png
% my model, response to deaths
S = zeros(T,1);
I = S;
R = S;
D = S;
C = S;
betat = S;
S(1) = N-1;
I(1) = 1;
I(1) = 1;
betat(1) = beta0;
for t = 1:T-1;
S(t+1) = S(t) -betat(t)*S(t)*I(t)/N;
I(t+1) = I(t) + betat(t)*S(t)*I(t)/N - gam*I(t);
R(t+1) = R(t) + gam*I(t)-thet*R(t);
D(t+1) = D(t) + delt*thet*R(t);
C(t+1) = C(t) + (1-delt)*thet*R(t);
betat(t+1) = exp(log(beta0) - alphaD*((D(t+1)-D(t))/N));
end;
figure;
%yyaxis left
plot((2:T)',[ I(2:T) 100*(D(2:T)-D(1:T-1)) ],'linewidth',2)
xlabel('Days')
ylabel('People');
axis([0 100 0 inf]);
yyaxis right
plot((1:T)',betat/gam,'--k','linewidth',2);
ylabel('Reproduction rate R0','color','k')
axis([0 100 0 inf]);
legend('Infected','100 x Deaths/day','R0 (right scale)','location','best')
title('BSIR model, R0 varies with death rate');
print -dpng std_sir_aD.png
