key: cord-0991226-zm728wk6
authors: Fernández-Villaverde, Jesús; Jones, Charles I.
title: Estimating and Simulating a SIRD Model of COVID-19 for Many Countries, States, and Cities
date: 2022-01-29
journal: J Econ Dyn Control
DOI: 10.1016/j.jedc.2022.104318
sha: ee9c033b1ac0e30f9d617955f4a7831d1ee72ca3
doc_id: 991226
cord_uid: zm728wk6

We use data on deaths in New York City, Madrid, Stockholm, and other world cities as well as in various U.S. states and other regions and countries to estimate, quickly and with limited data, a standard epidemiological model of COVID-19. We allow for a time-varying contact rate in order to capture behavioral and policy-induced changes associated with social distancing. We simulate the model forward to consider possible scenarios for various countries, states, and cities, including the potential impact of herd immunity on re-opening.

The sudden arrival of COVID-19 in the winter of 2020 highlighted the importance of estimating a standard epidemiological model of the epidemic quickly and with limited data. In this paper, we show how to tackle this challenge. We use data on deaths in New York City, Madrid, Stockholm, and other world cities as well as in various U.S. states, countries, and regions around the world during the first half of 2020 to estimate a SIRD model of COVID-19. Relative to existing frameworks, our contributions are:

• We do not use data on cases or tests because of differential selection in testing in different cities, states, and countries. Instead we only use data on deaths.

• We invert a standard SIRD epidemiological model and use the daily death series to recover a time-varying the basic reproduction number (i.e., the expected number of infections generated by one infection when all individuals are susceptible to infection) R 0t ≡ β t /γ to capture changes in behavior and policy that occur at different times and with different intensities in different locations. In essence, we apply a Solow residual approach: we assume the model fits the data exactly and back out the implied values of β t that make it so.

• We show how simulating our model after a location has reached a peak in the number of daily deaths results in very stable results going forward in time. In contrast, simulations of the future before a location reaches its peak are extremely noisy and sensitive to daily shocks.

• For simulations of future outcomes, we allow for feedback from daily deaths, d t , to future behavior according to R 0t = Constant·e −αdt as suggested by Cochrane (2020) . We estimate α from data for each country. There is tremendous heterogeneity across countries, so this parameter is not well-identified in our data. We estimate an average value of about α = 0.05 so that R 0 changes by 5% when daily deaths change by one and use this value in simulations of future outcomes.

• Our models allow us to back out the percentage of people who were infectious at the end of our sample as well as those who were ever infected versus those still susceptible; therefore, we can estimate the extent to which herd immunity effects are large. Given the epidemiological situation in mid-May 2020, we find moderate effects in New York City, noticeable effects in Italy, Sweden, and Spain, and negligible effects in New York state outside of New York City and in places like California.

We study a standard model of COVID-19 using common tools in econometrics, and then we analyze its main quantitative implications in ways that resemble how economists study other dynamic models. Our exercise can help us understand where a simple SIRD model has difficulties fitting observed patterns in the data and points out avenues for improvement while maintaining the virtues of simplicity and parsimony.

In the interest of space, we will report a very short summary of our results, up to mid-May 2020. By the end of May, the first wave of the epidemic was over in many cities, regions, and countries. Later waves of the epidemic need, to be analyzed in more detail, using models with time-varying parameters, such as the one in Arias et al. (2021) , and, consequently, much more powerful econometric techniques. Nonetheless, we have an online dashboard, https:// web.stanford.edu/~chadj/Covid/Dashboard.html, that reports data extended until October 9, 2020 for around 100 cities, states, and countries.

Much of the mathematical study of the spread of infectious diseases starts from the classic compartmental models of Kermack and McKendrick (1927) and Kermack and McKendrick (1932) . These models divide the population into several different compartments (e.g., susceptible, infective, recovered, deceased, ...) and specify how agents move across the separate compartments over time. The SIRD epidemic model that we analyze in this paper is one of the simplest of these compartmental models. Hethcote (2000) presents a useful overview of this class of models and some of their theoretical properties and Morton and Wickwire (1974) show how to apply optimal control methods to them.

The acute economic impact of the COVID-19 pandemic has generated a gigantic literature that we cannot review here except for pointing out a few papers that have particularly influenced our thinking (see Stock, 2020, and Avery et al., 2020 , for two general surveys of how economists have addressed this topic).

First, economists have argued that many of the parameters controlling the move among compartments are not structural in the sense of Hurwicz (1962) , but depend, instead, on individual decisions and policies. For example, the rate of contact that determines the number of new infections is a function of the endogenous labor supply and consumption choices of individuals. Hence, the rate of contact is amenable to being studied with standard decision theory models. See, for instance, Eichenbaum et al. (2021) and Farboodi et al. (2021) . Also, the recovery and death rates are not just clinical parameters, but can be functions of policy decisions such as expanded hospital capacity or priorities regarding the allocation of scarce ICU resources. Similarly, the case fatality ratio, a key figure to assessing the severity of the epidemic, is a complex function of clinical factors (e.g., the severity of a virus) and demographic and selection-into-disease mechanisms, which are themselves partly the product of endogenous choices (Korolev, 2020) . 1 Our paper builds on these ideas by allowing the infection rates to be influenced by social distancing and by letting many parameters vary across countries, states, and cities, which can proxy for demographic and policy heterogeneity.

Second, economists have been concerned with the identification problems of compartmental models. Many of these models are unidentified or weakly identified, with many sets of parameters that fit the observed data so far equally well but have considerably different long-run consequences. Atkeson (2020) and Korolev (2021) document this argument more carefully. Our findings corroborate this result and highlight the need to develop alternative econometric approaches.

Third, some researchers have dropped the use of compartmental models completely. Instead, they have relied on time-series models from the econometric tradition. See, for instance, Li and Linton (2021) and Liu et al. (2021) .

Let us close this section by pointing out that economists are pushing the study of compartmental models in a multitude of dimensions. Acemoglu et al. (2021) , Alvarez et al. (2021) , and Chari et al. (2020) characterize the optimal lockdown policy for a planner who wants to control the fatalities of a pandemic while minimizing the output costs of the lockdown. Berger et al. (2020) analyze the role of testing and case-dependent quarantines. Bodenstein et al. (2021) combine a compartmental model with a multisector dynamic general equilibrium model to capture key characteristics of the U.S. Input-Output Tables. Garriga et al. (2021) , Hornstein (2020) , and Toda (2020) study a variety of containment policies. More papers are appearing every day.

We follow standard notation in the literature. There is a constant population of N people, each of whom may be in one of five states:

The states -in temporal order-are

A susceptible person contracts the disease by coming into "adequate" contact with an infectious person, assumed to occur at rate β t I t /N , where β t is a time-varying contact rate parameter.

explicitly "causal" concept.

The starting value of β t , β 0 , reflects how the infection would progress if individuals behaved as they did before any news of the disease had arrived. We think of β 0 as capturing characteristics of the disease, fixed attributes of the region such as density, and basic customs in the region.

Over time, β t varies depending on how strong are the social distancing and hygienic practices that different locations adopt, either because of policy or simply because of voluntary changes in individual behavior. We will explain below how we recover β t from the data but, at this moment, we are not imposing any structure on its evolution. 2 The total number of new infections at a point in time is β t I t /N · S t . Infectiousness resolves at Poisson rate γ, so the average number of days a person is infectious is 1/γ: e.g., if γ = 0.2, a person is infectious on average for 5 days.

After the infectious period is over, a person is in the "Resolving" state, R. A constant fraction, θ, of people exit this state each period, and the case is resolved in one of two ways:

Death: fraction δ, Recovery:

fraction 1 − δ.

In preliminary work, we found it important to have a model that distinguishes between the infectious and the recovering periods. This distinction was key to matching the data with biologically plausible parameter values when we were putting restrictions on the time path of β t . It appears that the infectious period lasts on average about 4 to 5 days while cases take a total of about 2 to 3 weeks or even longer to resolve ( Bar-On et al., 2020) . 3 If one assumes people are infectious for this entire period, the model has trouble fitting the data.

The laws of motion related to the virus are then given by

We assume the initial stocks of deaths are set equal to zero. The initial stocks of infections and resolving cases, I(0) and R(0), are parameters that we will estimate.

Here we review the basic properties of this model when β t = β and the difference equations are replaced by differential equations (Hethcote, 2000) . A convention in epidemiological modeling is to recycle notation and let R 0 denote the basic reproduction number, that is, the expected number of infections generated by the first ill person when s 0 ≡ S 0 /N ≈ 1: More generally, if R 0 s 0 > 1, the disease spreads; otherwise, it declines quickly. One can see from this simple equation why R 0 > 1 is so natural: if people are infectious for 5 days and have lengthy contacts with even just two new people per day, for example, then R 0 = 10.

The initial exponential growth rate of infections is β − γ = γ (R 0 − 1). Another useful result concerns the long-run number of people who ever get infected (and therefore the fraction δ of these gives the long-run death rate). As t → ∞, the total fraction of people ever infected, e * , solves (assuming s 0 ≈ 1)

In other words, with a constant β, the long-run number of people ever infected is pinned down by R 0 ; the parameters γ and θ only affect the timing, holding R 0 constant. The long-run death rate is then δe * , which also depends only on R 0 (and δ). This explains why modeling the changing β associated with social distancing and better hygienic practices is so important. With a constant β, the initial explosion rate of the disease implies a value for β and then all the variables in the differential system are determined at that point. Instead, a changing β permits the initial exponential growth rate of deaths to be different from the long-run properties of the system, which is the point of adopting behavioral changes in society.

It turns out that recovering β t , a latent variable, from the data is straightforward without resorting to any complex filtering device.

We adopt the following timing convention. D t+1 is the stock of people who have died as of the end of date t + 1, so that ∆D t+1 ≡ d t+1 is the number of people who died on date t + 1 (daily deaths, in our estimating exercise).

We begin by using equation (4) to solve for various series involving R t+1 and its differences in terms of daily deaths:

Next, we use (3) and the expressions we just derived for R t+1 to solve for I t and its differences:

and applying the difference operator gives:

where ∆∆d t+3 ≡ ∆d t+3 − ∆d t+2 .

Taking the ratio of (9) to (8) gives:

Now, we can go back to our original SIRD model in equation (2) and rewrite it as

Solve this equation for β t by using equation (10) above to get:

.

This is one of the key equations in recovering β t . Notice, however, that this equation depends on S t . But since we have an initial condition for S 0 , we can use the SIRD model to get the updating equation for ∆S t+1 and we will be done. From (1) and using (8) to substitute I t :

Now, we only need to collect the last two equations together:

and:

With these two equations, an observed time series for daily deaths, d t , and an initial condition S 0 /N ≈ 1, we iterate forward in time and recover β t and S t+1 . Basically, we are using future deaths over the subsequent 3 days to tell us about β t today. While this means our estimates will be 3 days late (if we have death data for 30 days, we can only solve for β for the first 27 days), we can still generate an informative estimate of β t .

We can perform many exercises with the recovered β t . We can, for instance, simulate the model forward using the most recent value of β T and gauge where a region is headed in terms of the infection. And we can correlate the β t with other observables to evaluate the effectiveness of certain government policies such as mandated lockdowns.

Note, also, that β t determines the basic reproduction number, R 0t = β t × 1/γ under the prevailing social distancing and hygienic practices. We should be careful to distinguish this basic reproduction number from the effective reproduction number (i.e., the average number of new infections caused by a single infected individual at time t), which we will denote by R et . The latter considers the fraction of the population that is still susceptible. Since:

our procedure can also recover the effective reproduction number. This finding is interesting because this effective reproduction number is often reported by researchers due to the ease with which it can be estimated with standard statistical packages such as EpiEstim in R.

Now, we take our model to the data. The following parameters are assumed to be primarily biological and, therefore, fixed over time and the same in all countries and regions:

• γ = 0.2: In the continuous-time version of this model, the average length of time a person is infectious is 1/γ, so 5 days in our baseline. This choice is consistent with the evidence in Bar-On et al. (2020) . We also consider γ = 0.15 (7 day duration). The γ = 0.2 fit slightly better in our earlier work with more restrictions on β t , but it was not particularly well identified. 4

• θ = 0.1: In the continuous-time version of this model, the average length of time it takes for a case to resolve, after the infectious period ends, is 1/θ. With θ = 0.1, this period averages 10 days. Combined with the 5-day infectious period, this implies that the average case takes a total of 15 days to resolve. The implied exponential distribution includes a long tail that can be thought of as capturing the fact that some cases take longer to resolve.

• α = 0.05: For simulations of future outcomes, we allow for feedback from daily deaths per million people, d t , to future behavior according to R 0t = Constant · e −αdt as suggested by Cochrane (2020) . We estimate α i from data for each location i. There is tremendous heterogeneity across locations in these estimates, so a common value is not well-identified in our data. We estimate an average value of about α = 0.05 so that R 0t changes by 5% when daily deaths change by one. This is the value we use in simulations of future outcomes. More specifically, the mean value ofα i in location-specific regressions is 0.066 and the median value is 0.045. However, the standard deviation ofα i across locations is a very high 0.15. We report results with both α = 0 -i.e., assuming no feedback so that the final value of R 0t that we estimate in the data is assumed to hold in the future-as well as with α = 0.05. The presence of feedback is very clear in our estimation and strikes us as helpful to incorporate, so our baseline results below assume α = 0.05.

• δ = 1.0%: This parameter is crucial, and it would be great to have a precise estimate of it. Case fatality rates are not helpful, as we do not have a good measure of how many people are infected. Random testing for antibodies to detect how many people have ever been infected is quite informative about this parameter. We explain below how we use such data.

Seroepidemiological surveys. The most comprehensive evidence from the early stages of the COVID-19 epidemic we are aware of comes from a seroepidemiological national survey undertaken by the Spanish government from April 27 to May 11, 2020, to measure the incidence of SARs-CoV-2 in Spain. The survey was large, with 60,983 valid responses from individuals stratified in two stages. Combining the results from this survey with the measured sensitivity and specificity of the test, we conclude that the mortality rate of SARs-CoV-2 in Spain was between 1% and 1.1%. Because many of the early deaths in the epidemic were linked with mismanagement of care at nursing homes in Madrid and Barcelona that could have been avoided, we pick 1% as our benchmark value.

Since mortality rates are affected by the demographic composition of the population (with COVID-19 mortality rates increasing sharply with age), we obtained data on age distributions across countries from the U.N. population division. We decomposed the Spanish mortality rate by age, given the age-specific measured incidence of infection rates, and applied those age-mortality rates to the population shares of each country. To control for differences in life expectancy (and, hence, for the possibility that the age-specific mortality rate of an 80-year-old individual in a high life-expectancy country is equivalent to the age-specific mortality rate of a 70-year-old individual in a low life-expectancy country), we applied a correction based on the ratio of the life expectancy of each country with respect to Spain's life expectancy.

We found that, for most of the countries in our sample, the estimated mortality rate clusters around 1% (with or without the correction for life expectancy). For example, for the U.S., we found a death rate of 0.76% without correcting for life expectancy and 1.05% correcting for it. Therefore, and parsimoniously, we selected 1% as our baseline parameter value.

Other studies suggest similar values of δ. For instance, on April 23, 2020, Governor Andrew Cuomo announced preliminary results suggesting that 21% of New York City residents randomly tested from supermarkets and big-box stores had antibodies for COVID-19. According to the New York Department of Health (2020), it takes 3-4 weeks for these antibodies to form, so this suggests that around April 1, 21% of NYC residents were "ever infected." This infection rate is consistent with back-of-the-envelope calculations of death rates of around 0.8%-1.2%. Thus, we will report robustness results using death rates of 0.8% and 1.2%. 5

Data. Our data are taken from the GitHub repository of Johns Hopkins University CSSE (2020), which reports cumulative death numbers daily for countries, states, counties, and provinces throughout the world. The exception is for the international cities/regions of Lombardy, London, Madrid, Stockholm, and Paris. We obtain data for these locations from the various national vital statistics agencies.

Our sample for the results reported in this paper goes from the start of the epidemic until May 19, 2020. By the end of May, the first wave of the epidemic was over in many countries. When deaths are close to zero, our procedure often delivers negative values of R 0t : for example, small random changes from 1 death to 2 deaths a day imply second and third differences of the daily deaths that the standard SIRD model cannot rationalize. Also, our simple approach would need to be enriched to account for the repeated waves that arrived after the summer of 2020. This was done in follow-up research by one of us in Arias et al. (2021) . Nonetheless, our online dashboard reports data extended until October 9, 2020.

We manipulate the data in three important ways before feeding them into the model. First, on April 15, 2020, New York City added more than 3,500 deaths to its counts, increasing the total by more than 43%. We apply this same factor of proportionality (1.4325) to the deaths before April 15, 2020, to get a consistent time series for New York City. Second, The Economist (2020) reports that similar adjustments need to be made in other countries. In particular, vital statistics records in countries including Spain, Italy, England, France, and Sweden suggest that "excess deaths" relative to an average over past years exceed deaths officially attributed to COVID-19 by a large margin. Hence, we increase deaths in all non-New York City locations by 33% for all dates. 6 Finally, there are pronounced "weekend effects" in the raw data: there are days, often on the weekend or on a holiday, in the middle of the pandemic when a country reports zero deaths, only to make up for this with a spike in deaths in subsequent days. We initially ran the model with the raw data, and the model works fine. However, applying a 5-day centered moving average to the data produces more stable results, so we make this final adjustment.

Guide to Graphs. In the interest of space, we only report a small subset of our results. We invite the reader to check our detailed results on our online dashboard. In general, we will report cumulative deaths through the latest date, daily deaths (data and simulating forward), and cumulative deaths simulating forward. Data are shown as circles or bars, and simulations are solid lines. Each graph may have several lines, typically for one of two reasons. In some graphs, we show the simulations adding data from the last 7 days of our sample. This provides an intuitive assessment of how sensitive the simulations are to one or two recent observations. In other graphs, we show alternatives for baseline, "high," and "low" values of certain parameters. Figure 1 shows the estimates of R 0t = β t /γ for New York City. For the baseline parameter values, the estimates suggest that New York City began with R 0 = 2.7, so that each infected person passed the disease to nearly three others at the start. This estimate agrees with other findings and it is particularly plausible for such a high-density metropolitan area as New York City. 7 Social distancing is estimated to have reduced this value to below 0.5 by mid-April. After that, R 0t seems to fluctuate around 1.0. Figure 1 shows the estimates of R 0t = β t /γ for New York City. For the baseline parameter values, the estimates suggest that New York City began with R 0 = 2.7, so that each infected person passed the disease to nearly three others at the start.

This estimate agrees with other recent findings and its particularly plausible for such a high-density metropolitan area as New York City. 8 Social distancing is estimated to have reduced this value to below 0.5 by mid-April. After that, R 0t seems to fluctuate around 1.0.

It is worth briefly reviewing the data that allows us to recover R 0t . As discussed in Figure 1 : New York City: Estimates of R 0t = β t /γ It is worth briefly reviewing the data that allow us to recover R 0t . As discussed in Section 4, we invert the SIRD model and use the death data to recover a time series for R 0t such that the model fits the death data exactly. This inversion reveals that R 0t can be recovered from the daily number of deaths (d t+1 ), the change in daily deaths (∆d t+2 ), and the change in the change in daily deaths (∆∆d t+3 ). ing noise reduces precision), the latest estimate of R 0t we have for each location corresponds to May 9, even though our death data runs through May 19: we lose 2 observations for the moving average, 3 observations for the double differencing, and then truncate by an additional 5 days to improve precision. 9

Our estimation also allows us to recover the fraction of the population that is estimated to be currently infectious at each date. These results are shown for New York City in Figure 5 . For our baseline parameter values, this fraction peaks around April 1 at 5.7% of the population. By May 9, it is estimated to have declined 9 These graph of the underlying data are available for every location in our dataset in the extended results on our dashboard. Figure 4 shows the double difference. It is these HP-filtered data that are used in the construction of R 0t in Figure 1 . Because the HP filter has problems at the end of the sample (e.g., there are fewer observations so noise becomes more important, and double differencing noise reduces precision), the latest estimate of R 0t we have for each location corresponds to May 9, 2020, even though our death data run through May 19, 2020: we lose 2 observations for the moving average, 3 observations for the double differencing, and then truncate by an additional 5 days to improve precision.

Our estimation also allows us to recover the fraction of the population that is estimated to be infectious at each date. These results are shown for New York City in Figure 5 . For our baseline parameter values, this fraction peaks around April 1, 2020, at 5.7% of the population. By May 9, 2020, it is estimated to have declined to only 0.43% of the population. Figure 6 shows the time path of R 0t for several locations. There is substantial heterogeneity in the starting values, but they all fall and cluster around 1.0 once the pandemic is underway. By the end of our sample, the values of R 0 for Atlanta and Stockholm are noticeably greater than 1.0. Figure 7 shows the time path of the percentage of the population that is currently infectious, I t /N , for several locations. The waves crest at different times for different locations, and the peak of infectiousness varies as well. Table 1 summarizes these and other results for a broader set of our locations. The full table, together with around 39 pages of graphs for each location, is reported on our dashboard. Now is a good time to make a couple of general remarks about our estimation. First, as the number of daily deaths declines at the end of a wave -say for Paris, Madrid, and Hubei in the table- Table 1 summarizes these and other results for a broader set of our locations. The full table, together with ∼25 pages of graphs for each location, are reported on our dashboard. Now is a good time to make a couple of general remarks about our estimation. First, as the number of daily deaths declines at the end of a wavesay for Paris, Madrid, and Hubei in the table -the estimation of R 0t can become difficult and dominated by noise. In the extreme, for example, once total deaths are constant, our procedure gives β t = 0/0. One sign of such problems is that "today's" value of R 0 can fall to equal 0.20 -this is a lower bound that we impose Figure 6 : Estimates of R 0t = β t /γ the estimation of R 0t can become difficult and dominated by noise. In the extreme, for example, once total deaths are constant, our procedure gives β t = 0/0. One sign of such problems is that "today's" value of R 0 can fall to equal 0.20 -this is a lower bound that we impose on the estimation. When a location hits this lower bound, our routine ignores subsequent days of results because the model yields inconsistent result (e.g., negative new infectious). The notation "today" in the table refers to the last day for which we have results. Typically it is May 9, 2020, but in some cases it is earlier. Next, we turn to some general comments about the results. First, notice that the initial values for R 0 range from around 1.5 or lower in places like Minnesota, California, Norway, and Mexico to high values of 2.5 or more in major cities throughout the world. Second, the fraction of the population that is infectious at the peak is greater than 2% in the hardest-hit areas, but only reaches a maximum of 5.7% in New York City. Third, the fraction that is infectious at the end of the sample is typically lower. It has fallen below 0.4% in New York City (plus), Lombardy, Madrid, Paris, and Detroit but is greater than 0.7% in places including New Jersey, Stockholm, Philadelphia, and Chicago. It is even lower -below 0.1%-in the SF Bay Area, Washington state, and Germany. Finally, there is enormous heterogeneity in cumulative deaths per million people ("Total (pm) Deaths" in the table), both at the end of the sample and in the forward simulation for 30 days in the future (t+30). Figures 8, 9 , and 10 show how the model fits the New York City data for three values of δ: 0.01, 0.008, and 0.012. The main lesson is that the model fits the data very well with each of these parameter values: our procedure just adjusts the number of infected people to account for the same observed deaths. For example, with δ = 1.0%, our model implies that this number for April 1, 2020, was 17%. This compares very well with the observation that -as of April 20, 2020about 21% of New York City residents tested positive for antibodies of COVID-19 (New York Department of Health, 2020). Because antibodies only appear 3 to 4 weeks after infection, these antibody tests really tell us what the ever-infected rate was 3 to 4 weeks earlier. The subtitle lines for these three figures also report the "%Infected" at different dates. These are the percent of people who are estimated to have ever been infected with the virus. For New York City, the numbers as of early May are 26 percent, and then in 30 days are estimated to equal 27 percent, with a slightly higher value at the end of our simulation (the third number). We return in Section 7 to the implications of these high infection rates for herd immunity and re-opening.

The next set of graphs show results for Spain together with robustness to different The supertitle lines for these three figures also report the "%Infected" at different dates. These are the percentage of people who are estimated to have ever been infected with the virus. For New York City, the numbers as of early May 2020 are 26% percent, and then in 30 days they are estimated to equal 27%, with a slightly higher value at the end of our simulation (the third number). We return in Section 7 to the implications of these high infection rates for herd immunity and re-opening. Our dashboard reports similar exercises for many other locations. 

When we simulate the model for many countries and regions, we find two results. First, once countries or regions reach the peak and deaths start to decline, the forecasts converge well. Second, however, before that happens, the forecasts are very noisy. This makes sense: we are trying to forecast 30 to 60 days into the future based on 3 to 4 weeks of data using a very nonlinear model. We illustrate these points with the next two figures, Figures 11 and 12 , which show results at the end of our sample for New York City, now broadly defined to include the surrounding counties of Nassau, Rockland, Suffolk, and Westchester (which we call "New York City (plus)" in the graphs). In each figure, we see seven lines of forecasted daily and cumulative deaths. Each line corresponds to the forecast using one more day of observations. In both figures, the more recent observations push the forecast down (i.e., the top lines use fewer observations) and lowers its variance from day to day. This convergence of the forecast reflects how the first wave of COVID-19 was winding down in New York by late May 2020. Recall the role of the α feedback parameter. In the baseline simulation results, we assume R 0t = Constant · e −αdt where α = 0.05. This implies that if daily deaths rise, people adjust their behavior to reduce contacts, which reduces R 0t . Conversely, if daily deaths fall, people are more likely to go out and interact, which raises R 0t .

A point that is important to appreciate is that aggregating up from the city or county to the state and to the national level can be misleading. SIRD is a nonlinear model, so the results at the state level are not the same as the average of the results at the county level.

This point is easy to illustrate using data from New York. We report results for several different geographic regions. "New York City (plus)" includes New York City plus the four surrounding counties of Nassau, Rockland, Suffolk, and Westchester, with a total population of about 12 million. New York state is self-explanatory and has a population of about 20 million.

And "New York excluding NYC" is the difference between these other two: New York state excluding the NYC (plus) area, with a population of about 8 million. Now compare the results for these three regions, shown in Figures 13 through 15 . The results in New York state as a whole are driven entirely by New York City. For example, imagine (counterfactually) that there were no deaths outside of New York City. In this hypothetical case, deaths per million for New York state would look exactly like deaths per million for New York City, except scaled down by a factor of 12/20. Because of the lower deaths per million, the model would behave slightly differently. And yet New York outside of New York City could look very different. In fact, as the deaths in New York City decline, a potential rise in deaths outside of New York City could cause the state death numbers to exhibit a flattening or even a second peak. Another version of this same kind of geographic aggregation bias seems likely to occur for the United States itself. To see this, imagine 50 states that sequentially pass through the peak of daily deaths. The U.S. national number can be driven by New York (City!) for the first several weeks, then by New Jersey and Michigan, and then by Massachusetts and Pennsylvania. The U.S. graph may show a rise and then a very flat profile of deaths that persists for a long time before declining, as new regions within the country suffer through their peaks sequentially.

An important question at any stage of a pandemic is when to re-open the economy. The estimation we have conducted has something helpful to contribute to this point. Table 2 reports the estimated fraction of the population that had ever been infected as of May 9, 2020, for different countries and regions. Numbers for three different values of δ are also reported, with the baseline case of δ = 1.0% in the center column. Two key things stand out in the table. First, consider the baseline. As we discussed above, we estimate that 26% of New York City had ever been infected by late May 2020.

In contrast, only 4% of people in New York state outside of New York City and only 2% of Californians have ever been infected. There is enormous heterogeneity in ever-infected rates. Where do these numbers come from? In our model, the fraction δ of those infected eventually die, with the timing determined by γ and θ, but essentially suggesting that deaths at time t reflect infections from 15 days earlier. With an assumed death rate of δ = 1.0%, for each death, there are approximately 100 other people who have been infected. The large differences in the number of deaths per million in New York versus California then translate into these differences in infection rates. Interestingly, rates in Norway and South Korea are similarly very low, while ever-infected rates in Italy, Spain, and France are estimated to be around 6 to 8%. The second point is that these numbers are -in an obvious way-very sensitive to the assumed value of δ. If you double the death rate, you (roughly) halve the ever-infected rate. If you halve the death rate, you (roughly) double the infected rate. And as we discuss in more detail next, in thinking about herd immunity and re-opening the economy, knowing the fraction ever-infected is crucial, at least under the important assumption that antibodies give rise to immunity for an extended period of time.

There is an important complementarity here. We would like the death rate to be low, not just because it means that fewer people die, but also because it means that lots of people will already have been infected. For example, if the true death rate is 5 in 1000 rather than 10 in 1000, it means that 51% of New Yorkers had already been infected and the herd immunity effects would be very strong. In this sense, the finding that only 21% of New York City was ever infected as of April 1, 2020, was doubly bad news: it pushes up the death rate and means we are far from herd immunity, even in the place with the largest number of infections.

As Atkeson (2020) , Stock (2020) , and others have emphasized, random testing would have been extremely helpful in identifying which of these cases was relevant. Moreover, the table suggests that it was much more important to test in New York City than in California. So few people were likely infected in California that it would have been very hard to distinguish statistically between the different death rates, whereas even a few thousand random tests would have been very informative in New York City. This is a crucial point to remember for future epidemics.

This brings us to the next reason why knowing the percentage ever infected would be so useful. The complement of this number is the percentage of the population that is still susceptible to the virus at any given moment in time. Call this fraction s(t) ≡ S(t)/N (or better might be S(t)/(N − D(t)) but D(t) is so low that it makes no difference). 8 Recall from the basic SIR model that the virus will die out as long as R 0 (t)s(t) < 1, that is, if R 0t ≡ β t /γ is smaller than 1/s(t). The term s(t) is herd immunity. The fewer people who are susceptible and the more people who are recovered and hence immune, the less our random interactions result in infections. In particular, we can relax social distancing -increase β t and R 0t -to the critical value such that R 0t s(t) is just below one. That would mean that infected people infect fewer than one person on average, so herd immunity keeps the virus from re-surging. Table 3 shows these calculations for one month from the end of our sample (t + 30) given the baseline estimates from the model. For example, from the middle column, it is estimated that at t + 30, 78% of New York City (plus surrounding counties) would have still been susceptible. This means we could relax social distancing to the point where R 0 would rise to 1/0.78 = 1.3. This compares to the estimate for New York City at the end of the sample of 0.4 and the initial estimate of 2.6. In other words, New York City could move 41% ([1.3-0.4]/[2.6-0.4]) of the way back to normal and see no resurgence of the virus.

The rest of the state of New York, in contrast, is estimated to still have had 93% of the population susceptible a month from the end of our sample. So outside of the city, New York needed to maintain its R 0 at 1.1 -also its level at the end of our sample-to keep the virus from spreading. New York City and the rest of New York state needed different policies if the fraction of the population that remains susceptible was as different as these estimates imply. 9 Places with values of R 0 < 1 could have relaxed somewhat and still have kept the virus in check. But the basic news from this table is that with a death rate of 1%, there was very little accumulated herd immunity and that our scope for relaxing social distancing was limited (as shown by the later waves of the epidemic).

Finally, note that the SIRD model has "momentum." Even if an area has reached the threshold R 0 (t)s(t) < 1, we will continue to accumulate infections and deaths before the epidemic dies out fully. The number of these "overshoot" infections and deaths will depend on the number of infectious individuals when we reach R 0 (t)s(t) < 1. This observation is not a minor point. In a conventional SIRD model where R 0 (t) gives you herd immunity at 60% of the population, if we reach s(t) = 0.4 too fast, we can end up with over 90% of the population ever infected, that is, with an extra 30% of infections over those required to achieve herd immunity.

This means that we want to reach the threshold R 0 (t)s(t) < 1 or stay around it with very few infectious individuals to minimize "overshoot" infections. While setting up and solving an optimal control problem of the COVID-19 epidemic in the tradition of Morton and Wickwire (1974) to get to such an objective is beyond the scope of our paper, our empirical results can help to calibrate re-opening scenarios such as those quantitatively explored in Baqaee et al. (2020) . 

Our paper has presented a fast procedure to estimate a SIRD model with limited data. This exercise is particularly useful at the start of an epidemic, when a fast policy response is required and we cannot wait for months to implement more sophisticated econometric methods such as those in Arias et al. (2021) .

Relative to the standard SIRD model in the literature, we include a time-varying β, and therefore a time-varying R 0 . We invert the SIRD model to back out the daily values of R 0t that fit the death data. We see this as important for capturing behavioral changes by individuals in response to the pandemic as well as policy changes related to social distancing. We also include an additional "recovering" state that is consistent with the medical evidence that cases seem to be infectious for four to five days while taking a total of several weeks or more to resolve. These changes better connect the model to the epidemiology of the virus and are important in improving the model's ability to fit the data. Finally, we follow Cochrane (2020) and include feedback between R 0t and daily deaths in modeling the future of the epidemic. We hope that our empirical estimates will prove useful to others in thinking about the possible path that COVID-19 may take at different locations and in analyzimg future epidemics.

A Multi-Risk SIR Model with Optimally Targeted Lockdown

A Simple Planning Problem for COVID-19 Lock-down, Testing, and Tracing

Bayesian Estimation of Epidemiological Models: Methods, Causality, and Policy Trade-Offs

How Deadly Is COVID-19? Understanding the Difficulties with Estimation of Its Fatality Rate

Policy Implications of Models of the Spread of Coronavirus: Perspectives and Opportunities for Economists

Policies for a Second Wave

SARS-CoV-2 (COVID-19) by the Numbers

Testing and Reopening in an SEIR Model

Social Distancing and Supply Disruptions in a Pandemic

The Hammer and the Scalpel: On the Economics of Indiscriminate versus Targeted Isolation Policies during Pandemics

Characterizing the Reproduction Number of Epidemics with Early Subexponential Growth Dynamics

An SIR Model with Behavior

The Macroeconomics of Epidemics

Internal and External Effects of Social Distancing in a Pandemic

Optimal Management of an Epidemic: An Application to COVID-19

The Mathematics of Infectious Diseases

Social Distancing, Quarantine, Contact Tracing, and Testing: Implications of an Augmented SEIR-Model

On the Structural Form of Interdependent Systems

2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository

Deaths Reach 6 Times the Normal Level, Far More Than Coronavirus Count Suggests

A Contribution to the Mathematical Theory of Epidemics, Part I

Contributions to the Mathematical Theory of Epidemics. II -The Problem of Endemicity

What Does the Case Fatality Ratio Really Measure?

Identification and Estimation of the SEIRD Epidemic Model for COVID-19

When Will the Covid-19 Pandemic Peak?

Panel Forecasts of Country-Level Covid-19 Infections

On the Optimal Control of a Deterministic Epidemic

The NYSDOH Wadsworth Center's Assay for SARS-CoV-2 IgG

High Contagiousness and Rapid Spread of Severe Acute Respiratory Syndrome Coronavirus 2

Data Gaps and the Policy Response to the Novel Coronavirus

Tracking COVID-19 Excess Deaths Across Countries

Susceptible-Infected-Recovered (SIR) Dynamics of COVID-19 and Economic Impact