Notation

We denote \(X_{t,d}, d = 0, \dots, D\) as the number of cases occurring at time \(t\) which appear in the dataset with a delay of \(d\). For example, a delay \(d = 0\) means that a case occurring on day \(t\) arrived in the dataset on day \(t\). We only consider cases reporting within a maximum delay \(D\). The number of cases reporting for time \(t\) with a delay of at most \(d\) can be written as:

\[ X_{t, \le d} = \sum_{i=0}^d X_{t,i} \tag{1} \]

A special case of this is the “final” number of reported cases at time \(t\), denoted by \[ X_t = X_{t, \le D} = \sum_{i=0}^D X_{t,i} \]

For delays \(d < D\) we define the notation

\[X_{t,>d} = \sum_{i = d+1} ^{D} X_{t,i}\]

representing the number of cases still missing after \(d\) delay. In the following we use uppercase letters (\(X_t\)) for random variables, lower case (\(x_t\)) for the corresponding observations, and hats (\(\hat{x}_t\)) for estimated/imputed values. We refer to \(X_t\) to describe a random variable, \(x_t\) for the corresponding observation, and \(\hat{x}_t\) for an estimated/imputed value. The matrix \(\mathbf{x}\) with entries \(x_{t,d}, t = 1, \dots, t^*, d = 1, \dots, D\) is referred to as the reporting matrix. For a description of the nowcasting terms being used in this document (e.g. reporting triangle) and their corresponding abbreviations in the package, please consult the nowcasting nomenclature vignette.

In the case where \(t^*\) corresponds to the present date, all entries with \(t+d>t^*\) have yet to be observed and are thus still missing. As the available entries at its bottom form a triangle, this incomplete reporting matrix is referred to as the reporting triangle.

Pre-processing of the reporting triangle

All of the following steps require that the reporting triangle only has non-negative entries. In practice this is not necessarily the case. For instance, if the reporting triangle has been computed from increments in subsequent data snapshots, occasional downward corrections due to data entrance issues can cause negative entries. We therefore apply a pre-processing step to re-distribute negative entries across neighbouring cells with positive entries.

Delay distribution estimation

Estimating the delay distribution from a reporting triangle

In the case where \(t^*\) is the present day, such that only a reporting triangle with missing entries is available, the estimator \(\hat{\pi}_d\) from (2) can only be evaluated after discarding all data from the last \(D-1\) time points. In order to use these partial observations, we use a different representation of the delay distribution via terms of the form

\[ \theta_d = \frac{\pi_d}{\pi_{\le d-1}} \tag{3} \]

for \(d = 1, \dots, D\). Here, in analogy to (1) we write

\[ \pi_{\le d-1} = \sum_{d'=0}^{d-1} \pi_d. \tag{4} \]

The \(\theta_d\) can be estimated via

\[ \hat{\theta}_d = \frac{\sum_{t=t^* - N + 1}^{t^*-d} x_{t, d}}{\sum_{t= t^* - N + 1}^{t^*-d} x_{t, \leq d - 1}}, \]

and translated to estimates \(\hat{\pi}_0,..., \hat{\pi}_D\) via the recursion

\[ \hat{\pi}_{\leq d} = (1+\hat{\theta}_d)\hat{\pi}_{\leq d-1} \]

subject to the constrain that \(\sum_{d = 0}^D \pi_d = 1\).

We note that this method is equivalent to the so-called chain ladder method^[2], adapted to our notation in terms of reporting triangles (rather than development triangles as used in accounting). The intuition behind the method relies on the basic assumption is that the ratio between the delays in the complete observations is equivalent to the ratio between the delays in the partially observed components of the matrix (Figure 1).

Figure 1: Visual description of the iterative “completing” of the reporting triangle, which relies on the assumption that the ratio between delays for fully observed reference times is consistent with the ratio between those same delays in the partially observed data.

The delay_estimate() function ingests either a reporting matrix or a reporting triangle and uses the last n rows to compute an empirical delay probability mass function (PMF), returning a simplex vector indexed starting at delay 0.

Point nowcast generation

We now address the computation of a point nowcast, i.e., expected final case numbers \(\hat{x}_t, t = t^* - D + 1, t^*\). These are based on the reporting triangle, more specifically the preliminary totals \(x_{t, \leq t^* - t}\), and the estimated delay distribution, \(\hat{\pi}_d, d = 0, \dots, D\). In the following we will denote the current time \(t^*\) as the nowcast time, and the time \(t = t^* - D, \dots, t^*\) as the reference time. The difference \(j = t^* - t\) will be called the horizon.

An intuitive approach, used in^[1] and the standard chain ladder technique, is to simply inflate the current total for a reference time \(t\) by the inverse of the respective probability of observation up to time \(t^*\),

\[ \hat{x}_t =\frac{x_{t, \leq t^*-t}}{\pi_{\leq t^*-t}}. \]

This, however, is not well-behaved if no cases at \(t\) have been observed yet, i.e., \(x_{t, \leq t^* - t} = 0\). Then \(\hat{x}_t\) is likewise zero, which yields problems in our uncertainty quantification method (see next section). Motivated by a Bayesian argument (see Zero-handling below) we therefore use the expression

\[ \hat{x}_t = \frac{x_{t, \leq t^*-t} + 1 - \pi_{\leq t^*-t}}{\hat{\pi}_{\leq t^*-t}} \tag{5} \]

instead. This yields essentially identical results for large \(x_{t, \leq t^* - t}\), but produces positive \(\hat{x}_t\) even for preliminary zero values \(x_{t, \leq t^* - t} = 0\).

For our uncertainty quantification scheme we require not only estimated totals \(\hat{x}_t\), but all entries \(\hat{x}_{t,d}\) of a point nowcast matrix. For \(t = t^* - D + 1, \dots, t^*, d > t^* - t\) these are obtained as

\[ \hat{x}_{t,d} = \hat{\pi}_d \times \hat{x}_t. \]

The apply_delay() function ingests a reporting triangle and a delay PMF and returns a point nowcast matrix by “filling in” the missing elements of the reporting triangle. The generation of a point nowcast is described in the schematic below (Figure 2), demonstrating how the missing elements of the reporting triangle are estimated and can then be used to generate a point estimate of the final counts at each reference time.