
figuresection
Derivations for the dynamic models for \(\lbrace R_t\rbrace \)
Here we define \(L_t=\log (R_t)\)  for simplicity of notation and neglect the first time points (before March 7 2020) for which \(R_t\)  is constant. Note that for each model, the sufficient statistics involved are all easily updated with increasing \(t\) .
The random walk model
Assume a prior \(p(\tau )=\text{Gamma}(\alpha _\tau ,\beta _\tau )\)  with \(\tau =1/\sigma _R^2\) .
Given \(\lbrace R_t\rbrace \) , the posterior for \((\phi _c,\tau )\)  becomes
\(p(\tau |\mathbf {R}_{1:t})\propto &\tau ^{\alpha _\tau -1}\exp (-\beta _\tau \tau ) \prod _{s=1}^t\tau ^{1/2}\exp (-\tfrac{1}{2}\tau (L_s-L_{s-1})^2)\\\propto &\tau ^{\alpha _\tau +0.5t}\exp (-[\beta _\tau +0.5\sum _{s=1}^t(L_s-L_{s-1})^2]\tau )\\\propto &\text{Gamma}(\tau ;\alpha _\tau +0.5t,\beta _\tau +0.5\sum _{s=1}^t(L_s-L_{s-1})^2).\) 

Simulation of \(L_{t+1}\)  can then be performed by first simulating
\(\tau \)  from the Gamma distribution and thereafter
\(R_{t+1}\)  from the model.

The piecewise constant model
Assume a prior \(p(\phi _c,\tau )=\text{Beta}(a_\phi ,b_\phi )\text{Gamma}(\alpha _\tau ,\beta _\tau )\)  with \(\tau =1/\sigma _R^2\) .
Given \(\lbrace R_t\rbrace \)  and changepoint indicators \(C_t\)  (\(=1\)  if there is a changepoint at time \(t\) ), the posterior for \((\phi _c,\tau )\)  becomes
\(p(\phi _c,\tau |\mathbf {R}_{1:t})\propto &\phi _c^{a_\phi -1}(1-\phi _c)^{b_\phi -1}\tau ^{\alpha _\tau -1}\exp (-\beta _\tau \tau )\times \\&\prod _{s=1}^t[(1-\phi _c)^{1-C_s}\phi _c^{C_s}\tau ^{C_s/2}\exp (-\tfrac{1}{2}\tau C_s(L_s-L_{s-1})^2)]\\\propto &\phi _c^{a_\phi +n^c_t-1}(1-\phi _c)^{b_\phi +t-n_t^c-1}(1-\phi _c)^{n_t^c}\times \\&\tau ^{\alpha _\tau +0.5n_t^c}\exp (-[\beta _\tau +0.5\sum _{s=1}^tC_s(L_s-L_{s-1})^2]\tau )\\\propto &\text{Beta}(\phi _c;a_\phi +n_t^c,b_\phi +t-n_t^c)\times \\&\text{Gamma}(\tau ;\alpha _\tau +0.5n_t^c,\beta _\tau +0.5\sum _{t=1}^TC_s(L_s-L_{s-1})^2)\) 

where \(n^c_t\)  is the number of changepoints up to timepoint \(t\) . Simulation of \(R_{t+1}\)  can then be performed by first simulating
\((\phi _c,\tau )\)  from the Beta-Gamma distribution and thereafter
\(R_{t+1}\)  from the model.

The autoregressive model
We consider here only the case where \(\mu =0\) .
Assume
\(p(a,\tau )=&N(a_0,[\kappa _0\tau ]^{-1})\times \text{Gamma}(\alpha _\tau ,\beta _\tau )\) 

Note first that
\(g_t(a)\equiv &\kappa _0(a-a_0)^2+\sum _{s=1}^t(L_s-aL_{s-1})^2\\=&\kappa _0 a^2-2\kappa _0 a_0a+\kappa _0a_0^2+\sum _{s=1}^tL_s^2-2a\sum _{s=1}^tL_sL_{s-1}+a^2\sum _{s=1}^tL_{s-1}^2\\=&(\kappa _0+\sum _{s=1}^tL_{s-1}^2)a^2-2[\kappa _0 a_0+\sum _{s=1}^tL_sL_{s-1}]a+\kappa _0 a_0^2+\sum _{s=1}^tL_s^2\\\ =&(\kappa _0+\sum _{s=1}^tL_{s-1}^2)(a-\hat{a}_t)^2-(\kappa _0+\sum _{s=1}^tL_{s-1}^2)\hat{a}_t^2+\kappa _0 a_0^2+\sum _{s=1}^tL_s^2\\\multicolumn{2}{l}{\text{with}}\\\hat{a}_t=&\frac{\kappa _0 a_0+\sum _{s=1}^tL_sL_{s-1}}{\kappa _0+\sum _{s=1}^tL_{s-1}^2}.\) 

This gives
\(p(a|\tau ,R_{1:t})\propto &\exp \left(-0.5\tau g_t(a)\right)\\\propto & N\left(\hat{a}_t,\tau ^{-1}[\kappa _0+\sum _{s=1}^tL_{s-1}^2]^{-1}\right).\) 

Further, for any value of \(a\) , we have
\(p(\tau |R_{1:t})\propto &\frac{p(\tau )p(R_{1:t}|\tau ,a)}{p(a|\tau ,R_{1:t})}\\\multicolumn{2}{l}{\text{and using $a=\hat{a}_t(a)$, we obtain}}\\p(\tau |R_{1:t})\propto &\frac{\tau ^{\alpha _\tau -1}\tau ^{0.5t}\exp (-\beta _\tau \tau )\exp (-0.5\tau g_t(\hat{a}))}{\tau ^{0.5}}\\\propto &\text{ Gamma}(\alpha _\tau +0.5t,\beta _\tau +0.5g_t(\hat{a}_t)).\) 


Parameter settings
The full model includes several parameter specifications that are fixed throughout all runs. Some sensitivity analysis is provided in Section .

\(R_t\)  was set to a constant \(R_0\)  from February 17 2020 to March 8 2020. We estimate \(R_0\)  and used a \(N(1.192,0.1)\)  prior for \(\log R_0\)  (giving an expectation of \(R_0\)  equal to 3.3, which is compatible with the literature).

The amplification factor for imported cases is set to 2.8 from [1]} .

In order to take into account regional variability, the number of hospitalised individuals is assumed to follow a beta-binomial distribution with expectations equal to the probabilities given in Table REF  and with a shape parameter equal to 8. The time-dynamic changes in the probabilities are based on [1]}.

<TABLE>

The distribution for the delay from infected to hospitalised is assumed to follow a negative binomial distribution with parameters \((8.87,3.40)\)  until August 1 2020 and thereafter with parameters \((7.60,3.34)\) , as estimated from registry data.

The distribution of the time between infection and testing positive is assumed to follow a discrete distribution with probabilities \(0.05,0.080,0.16,0.4,0.3,0.01\)  for delays of 1-6 days.

In the model (REF ) for the detection probability, we used \(\pi _0=-1.017,\pi _1=0.00013\) . The number of detected individuals was assumed to follow a beta-binomial distribution with expectations given by ( REF ) and a shape parameter equal to 8.


Additional results
Figure REF  explores the sensitivity of results with respect to some of the parameters given in Section . In this figure (where only posterior medians are shown), the black line corresponds to weekly averaged fixed lag estimates of \(\lbrace R_t\rbrace \)  based on the autoregressive model (the middle panel in Figure REF ). For the red curve, we changed the prior for \(\tau =1/\sigma ^2\)  to a \(\text{Gamma}(1.2,0.14)\)  distribution.
For the blue curve, we changed the values of \((\pi _0,\pi _1)\)  to (-1.175,0.000074) (estimates obtained from another model on the same data). The green curve corresponds to a change in the shape parameters in the negative binomial distribution for hospitalisation and test data from 8 to 3. In all cases the results are stable, which was also the case when we compared 50% credibility intervals (not shown).
Figure REF  shows the fixed-24-days-lag smoothed estimates of daily \(R_t\)  values for the three models. This plot is not averaged over each preceding week, as is done in Figure REF  in the main text.
Figures REF -REF  show (filtered) parameter estimates of the dynamics of \(\lbrace R_t\rbrace \)  for the three models considered. All parameter estimates are based on the use of sufficient statistics [3]}, [4]} in the SMC algorithm.
Table REF  shows that although there is some Monte Carlo variability in the estimation of the marginal likelihoods, the differences between the models are statistically significant.

Figure REF  summarises results from the simulation experiments reported in Section REF . For each of the 20 simulation experiments (with the same underlying true \(\lbrace R_t\rbrace \)  process), posterior medians were obtained similar as to the runs on the real data. The figure shows the variability in these medians over different data sets, demonstrating the strong information content in such data about the \(\lbrace R_t\rbrace \)  process.
Figure REF  gives a comparison of one run based on \(B=20\,000\)  particles and one run with \(B=2\,000\)  particles, indicating that the number of particles used are sufficient for reliable estimates of the \(\lbrace R_\rbrace \)  process.
<TABLE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE>

Estimation of hyper-parameters
As described in [5]}, there are many different approaches to estimate static hyper-parameters, both online and offline.
In the paper we focused on online methods for estimation of parameters in the \(\lbrace R_t\rbrace \)  process, based on the use of sufficient statistics [3]}, [4]}. As pointed out by one of the referees, this method can suffer from degeneracy and converge towards different limit distributions in different runs. Such degeneracy problems can in particular occur for long time series. We validate the estimates obtained in two ways. We concentrate on the AR process, both due to the fact that this is our best model and that it turned out to be the most problematic one, with respect to parameter estimation.
Figure REF  shows the parameter estimates obtained from two independent runs. The agreement of these runs indicates that the estimation procedure performs quite well in our case. We do however emphasise that even though the problems with degeneracy did not seem to be present for the model considered here, if we considered an AR process where also \(\mu \)  is estimated, indeed different runs of the SMC algorithm resulted in somewhat different parameter estimates. The problems in that case occurred on the day where test data were introduced.
<FIGURE>An important feature of SMC methods is that unbiased estimates of marginal likelihoods are available, making the pseudo-marginal approach by [8]} possible. Given a current parameter \(\mathbf {\theta }\)  with corresponding estimate of the marginal likelihood \(\hat{p}(\mathbf {y}_{1:T};\mathbf {\theta })\) , a proposal \(\tilde{\mathbf {\theta }}\sim q(\tilde{\mathbf {\theta }}|\mathbf {\theta })\)  is generated and an estimate \(\hat{p}(\mathbf {y}_{1:T};\tilde{\mathbf {\theta }})\)  is obtained by a new run of the SMC algorithm. The new proposal will then be accepted with probability
\(\min \lbrace 1,\frac{p(\tilde{\mathbf {\theta }})\hat{p}(\mathbf {y}_{1:T};\tilde{\mathbf {\theta }})q(\mathbf {\theta }|\tilde{\mathbf {\theta }})}{p(\mathbf {\theta }\hat{)} p(\mathbf {y}_{1:T};\mathbf {\theta })q(\tilde{\mathbf {\theta }}|\mathbf {\theta })}\rbrace ,\) 
corresponding to the Particle Metropolis-Hastings Algorithm in [9]}.
Practical use of this algorithm requires tuning of the proposal distribution \(q(\tilde{\mathbf {\theta }}|\mathbf {\theta })\)  and of the number of particles used in each SMC run.
Concerning the latter, the number of particles in each SMC run needs to be reduced dramatically compared to a single run of the SMC algorithm since these runs have to be repeated many times. Although [9]} showed that the number of particles should increase with time, we have only considered a fixed number of particles in our testing. Even with many trials on the tuning parameters, we struggled to obtain convergence within reasonable time (running the algorithm for several days). However, we stored all the proposed parameter settings together with their likelihood values. We then fitted a two-dimensional GAM model to the log-likelihood values using \(a\)  and \(\sigma \)  as covariates. More than 500 samples of \((a,\sigma )\)  were obtained through multiple runs of the PMCMC algorithm with different starting points (running in total for several days). Figure REF  shows the fitted GAM curve with the red circles corresponding to the samples obtained from the SMC algorithm. The SMC samples seem to be within the higher likelihood region of the fitted curve. Note that due to lack of converge, the PMCMC algorithm has probably not been able to explore the full posterior properly. Further, the "observed" values displayed as small dots in the plot are all proposed values, not only the accepted ones, so these points should not be used to reflect the whole distribution.
<FIGURE>