This is an attempt to reproduce Figure 1 from Imai et al. 2020, using code from Riou and Althaus (code here). I have included a Jupyter notebook (which includes the output), as well as RMarkdown and plain R files, generated using Jupytext.
Details missing from Imai et al. that I had to guess:
- The initial time (I took it to be 2019-12-01).
- The distribution of generation times (I took it to be gamma with standard deviation of 3.8).
- The summary curve is generated using the median of all cases at a given time (similarly with the 95 percentile range). The language in the figure legend is ambiguous, and could mean generating a median trajectory.
I have also tweaked the Riou and Althaus code, changing the parameterization and wrapping the main loop into a function, and have also used a fixed seed. To reduce computational burden, I have only generated 500 trajectories (cf 5000 in Imai et al.), but this should not affect the results (particularly the median curve) too much. The details of the implementation are missing from Imai et al., and so they could well be using a different language, random number generator, etc..
At present, the figures do not match up; the median number of cases should go through 4000 cases on 2020-01-18. This discrepancy may well be due to my errors in the above assumptions and code refactoring.
So this is a bit faster; the use of
data.table
is not the main contributor to that. Also set it back to 5k samples.