Gompertz curve

Our model is based on the sigmoidal growth curve pattern that characterizes epidemics. We start by assuming that the cumulative number of cases follows a Gompertz curve, which is a parsimonious model for the canonical sigmoid shape of cumulative case numbers in a one-wave epidemic. Over the course of a wave, the number of new infected cases increases up to a peak before declining to zero as the pool of susceptible individuals declines. Specifically, if the cumulative number of cases at time \(t\), \(\mu(t)\), follows a Gompertz curve, we can write \[\mu(t) = \bar{\mu}\exp\{\gamma_0 e^{\gamma t}\},\] where \(\bar{\mu}\) is the unknown saturation level for the cumulative number of cases, \(\gamma_0 < 0\) is a parameter related to \(\mu(0)\) and \(\gamma < 0\) is the growth rate parameter. Defining \(\dot{\mu}(t) \equiv d\mu(t)/dt\) and \(g(t) \equiv {\dot\mu(t)}/{\mu(t)}\), it is straightforward to show that

\[\begin{align} \ln g(t) &= \delta + \gamma t, & (1) \end{align}\]where \(\delta = \ln \gamma_0 \gamma\).

The observational model needs to be specified in discrete, rather than continuous, time. This is straightforward. Let \(Y_t\) be the observed cumulative number of cases on day \(t\) and \(y_t = Y_t-Y_{t-1}\) be the number of daily new cases.¹ We can then define the growth rate of \(Y_t\) as \(g_t = y_t/Y_{t-1}\) and replace \(\ln g(t)\) with \(\ln y_{t}-\ln Y_{t-1}\).

Dynamic Gompertz model

The deterministic trend implied by (1) is too inflexible for practical time-series modeling of an epidemic. Replacing it with a stochastic trend allows the model to adapt to changes in dynamics during the course of the epidemic. We call this stochastic-trend counterpart of (1) the dynamic Gompertz model. It is a local linear trend model specified as

\[\begin{align} \ln g_{t}&=\delta_{t}+\varepsilon_{t}, \;\;\;\;\;\varepsilon_{t}\sim NID(0,\sigma_{\varepsilon }^{2}),\;t=2,...,T, & (2) \end{align}\] where \(\ln g_{t}=\ln y_{t}-\ln Y_{t-1}\) and

\[\begin{align} \delta_{t} &=\delta_{t-1}+\gamma_{t-1}, & (3)\\ \gamma_{t} &=\gamma_{t-1}+\zeta_{t},\;\;\; \zeta_{t}\sim NID(0,\sigma_{\zeta }^{2}), & (4) \end{align}\] where the disturbances \(\varepsilon_{t}\) and \(\zeta_{t}\) are mutually independent, and \(NID(0,\sigma^2)\) denotes normally and independently distributed with mean zero and variance \(\sigma^2\). Note that the larger the signal-to-noise ratio, \(q_{\zeta }=\sigma_{\zeta }^{2}/\sigma_{\varepsilon }^{2}\), the faster the estimate of the slope parameter, \(\gamma_t\), which can be interpreted as the growth rate of the growth rate of cumulative cases, changes in response to new observations. Conversely, a lower signal-to-noise ratio induces more smoothness to the estimates. When \(\sigma_{\zeta }^{2}=0\), the trend is deterministic as in (1).

State space form and estimation

It is convenient to write the dynamic Gompertz model in general state space form: \[\begin{align*} \ln g_t &= Z\alpha_t + \varepsilon_t & \varepsilon_{t} \sim NID(0,\sigma^2_{\varepsilon})\\ \alpha_{t+1} &= T\alpha_{t} + R \eta_t & \eta_{t} \sim NID(0,Q) \end{align*}\] with \[\begin{equation*} \alpha_t = \left(\delta_t, \gamma_t\right)', \; Z = \left(1, 0\right), \; \eta_t = \left(0, \zeta_t\right)', \; T = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, \; R = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}, \; Q = \begin{pmatrix} 0 & 0 \\ 0 & \sigma^2_{\zeta} \end{pmatrix}. \end{equation*}\] This model can be estimated using techniques based on the Kalman filter once a prior is specified. The prior is

\[\begin{equation*} (\delta_1, \gamma_1)' \sim N(a_1,P_1), \end{equation*}\]

where \(a_1\) is a \(2\times 1\) vector of prior means and \(P_1\) a \(2\times 2\) prior variance matrix. We use a diffuse prior due to the absence of prior information about the epidemic when the model is first estimated: i.e., we set \(a_1 = (0,0)'\), \(P_1=\kappa I\), and let \(\kappa \to \infty.\) Model estimation, including implementation of the diffuse prior, is carried out using the KFAS package Helske (2017).

The Kalman filter outputs estimates of the state vector \((\delta_{t},\gamma_{t})^{\prime }.\) The estimates at time \(t\) conditional on information up to and including time \(t\) are denoted \((\hat\delta_{t \mid t},\hat\gamma_{t \mid t})^{\prime }\) and given by the contemporaneous filter; the predictive filter estimates the state at time \(t+1\) from the same information set, outputting \((\hat\delta_{t+1 \mid t},\hat\gamma_{t+1 \mid t})^{\prime }.\)

It may be useful to review past movements of the state vector \((\delta_{t},\gamma_{t})^{\prime }.\) This can be done using the smoothed estimates \((\hat\delta_t,\hat\gamma_t)^{\prime}\), which denotes the estimates of the state vector at time \(t\) based on all \(T\) observations in the series.

Estimation of the unknown variance parameters (\(\sigma^2_{\varepsilon}\) and \(\sigma^2_{\zeta}\)) is by maximum likelihood (ML) and is carried out using KFAS following the procedure described in Helske (2017). We retain the option of either estimating the signal-to-noise ratio \(q_{\zeta}\), or of fixing it at a plausible value. In practice, for coronavirus applications, we set the value of \(q_{\zeta}\) based on experience and judgment, reducing the number of parameters to be estimated by one. Tests for normality and residual serial correlation are based on the standardized innovations, that is one-step ahead prediction errors, \(v_{t}=\ln g_{t}-\delta_{t \mid t-1},\) \(t=3,...,T.\)

Daily effects, which are generally quite pronounced in the coronavirus data, can be included in the model as described in the Appendix.

Forecasts and peak prediction

Forecasts of future observations are obtained from the predictive recursions

\[\begin{align} \nonumber \widehat{g}_{T+\ell \mid T} &=\exp (\hat\delta_{T\mid T}+\hat\gamma _{T\mid T}\ell ),\;\;\;\;\ell =1,2,...\\ \widehat{\mu }_{T+\ell \mid T} &=\widehat{\mu }_{T+\ell -1\mid T}(1+ \widehat{g}_{T+\ell \mid T}) \end{align}\] so that \[\begin{align} \widehat{y}_{T+\ell \mid T}&=\widehat{g}_{T+\ell \mid T}\widehat{\mu } _{T+\ell -1\mid T} = Y_T \exp\hat\delta_{T+\ell|T} \prod_{j=1}^{\ell-1}(1+\exp \hat\delta_{T+j|T}) & (5) \end{align}\] and \(\widehat{Y}_{T+\ell \mid T}=\widehat{\mu }_{T+\ell \mid T};\) the initial value is \(\widehat{\mu }_{T\mid T}=Y_{T}.\)

We construct state-implied forecast bands for \(y_t\) by propagating the predictive distribution of the latent state vector. In the one-step-ahead case the implied distribution for \(y_{T+\ell}\) is lognormal, but beyond one step the cumulative structure of the model means that the forecast path is path-dependent. The bands should therefore be interpreted as state-implied forecast bands rather than as calibrated prediction or confidence intervals.²

The filtered growth rate \(\hat{g}_{y,t\mid t}\) of new cases \(y_{t}\), can be extracted from the continuous-time incidence curve: \(\mu^{\prime}(t) = g(t) \mu(t)\), where \(\mu(t)\) is the growth curve and \(g(t)\) is its growth rate. Taking logarithms and differentiating we get

\[\begin{align} \hat{g}_{y,t\mid t}&=\hat{g}_{t\mid t}+\hat\gamma_{t\mid t}, & (6) \end{align}\]

where \(\hat{g}_{t\mid t}=\exp \hat\delta_{t\mid t}.\) The sampling variability of \(\hat{g}_{t\mid t}\) is dominated by that of \(\hat\gamma_{t\mid t}\) (see Harvey and Kattuman (2021)). Therefore when constructing confidence intervals for \(\hat{g}_{t\mid t}\) we treat \(\hat{g}_{y,t}\) as if it has a normal distribution centered on \(\hat{g}_{y,t\mid t}\) with variance \(\text{Var}(\hat\gamma_{t\mid t})\).

Even when the nowcast \(\hat{g}_{y,T\mid T}\) is positive and daily cases are growing, there will be a saturation level for the cumulative total, \(Y_{t},\) so long as \(\hat\gamma_{T\mid T}\) is negative. The nowcasts of \(y_{t}\) peak when \(\hat{g}_{y,t\mid t}=0,\) which requires \(\hat\gamma_{t\mid t}\) to be sufficiently negative to outweigh \(\hat{g}_{t\mid t},\) which is, of course, always positive. This can be seen from the expression for the growth rate of daily cases: \[\begin{align} \hat{g}_{y,T\mid T}&=\exp \hat\delta_{T\mid T}+\hat\gamma_{T\mid T}=\hat{g}_{T\mid T}+\hat\gamma_{T\mid T}. & (7) \end{align}\]

When \(\hat\gamma_{T\mid T}\) is negative, there is a flattening of the curve and a signaling of an upcoming peak in the trend of \(y_{t}.\) As shown in [HK, p10], the peak in the trend is predicted to be \(\ell _ {T}\) days ahead where³

\[\ell_{T}=\frac{\ln (-\hat\gamma_{T\mid T})-\hat\delta_{T\mid T}}{\hat\gamma_{T\mid T}}=\frac{\ln (-\hat\gamma_{T\mid T}/\hat{g}_{T\mid T})}{\hat\gamma_{T\mid T}},\;\;\;\; -\hat{g}_{T\mid T}<\hat\gamma_{T\mid T}<0. \]

The generation of forecasts is demonstrated in Section 4.

Reproduction Number \(R_t\)

The path of the epidemic is best tracked by nowcasts and forecasts of \(g_{y,t}\), the growth rate of \(y_t\), which are constructed by HK from the filtered estimates in the state space model, (2), (3) and (4). Wallinga and Lipsitch (2007) describe how the estimates of \(g_{y,t}\) can be translated into estimates of the instantaneous reproduction number \(R_{t}.\) Harvey and Kattuman (2021) propose \[\begin{align} \widetilde{R}_{t,\tau}&=1+\tau g_{y,t\mid t}\;\;\;\text{or}\;\;\;\widetilde{R}_{\tau,t}^{e} =\exp(\tau g_{y,t\mid t}),\;\;\; & (8) \end{align}\] where \(\tau\) is the generation interval – the typical number of days between an infected person becoming infected and them transmitting the disease to someone else. We construct approximate pointwise uncertainty bands for \(\widetilde{R}_{t,\tau}\) and \(\widetilde{R}_{\tau,t}^{e}\) by substituting the upper and lower bounds of the pointwise uncertainty bands for \(g_{y,t}\) into (8). These are approximate bands for the latent reproduction-number summaries implied by the fitted state-space model, conditional on the chosen generation interval \(\tau\); they are not Bayesian credible intervals. See Harvey et al. (2021) for an application.

The estimates of \(R_{t}\) can be used for tracking and forecasting the epidemic. The nowcasts of \(y_{t}\) peak when \(\hat{g}_{y,t\mid t}=0\), corresponding to \(\widetilde{R}_{t,\tau}=\widetilde{R}_{\tau,t}^{e}=1.\) Based on (7), predictions of \(g_{y,t}\) are given by

\[\begin{align} \hat{g}_{y,T+\ell\mid T} & =\exp\hat\delta_{T+\ell\mid T}+\hat\gamma_{T+\ell\mid T}=\exp(\hat\delta_{T\mid T}+\hat\gamma_{T\mid T}\ell)+\hat\gamma_{T\mid T},\;\;\;\;\;\ell=1,2,... & (9) \end{align}\]

We can then obtain predictions of \(R_{t},\) as in (8). If \(\hat\gamma_{T\mid T}\) is zero, the estimated growth of \(y_{t}\) is exponential and it is helpful to characterize it by the doubling time, \(\ln2/\hat{g}_{y,T\mid T} =0.693\exp(-\hat\delta_{T\mid T}).\)

When \(\exp \hat\delta _{T\mid T}+\hat\gamma _{T\mid T}>0\), the nowcast \(\hat{g}_{y,T\mid T}\) is positive and the estimate of \(R_{t}\) given by (8) is greater than one. So long as \(\hat\gamma _{T\mid T}\) is negative, then as \(T\rightarrow \infty ,\) \(\widetilde{R}_{\tau ,T+\ell \mid T}^{e}\rightarrow \exp (\tau \hat\gamma _{T\mid T})<1,\) and a saturation level for \(Y\) appears on the horizon.

We now turn to case where \(\gamma _{t}\) potentially turns positive in a typically short-lived phase, as a new wave emerges.

Reinitializing the data series

In practice the reinitialization rule is implemented in two stages, following a real-time trigger condition. Let \(T\) denote the current forecast origin and let \(\mathcal{F}_T\) denote the information available up to and including date \(T\). Write \(P_{\gamma,T\mid T} = Var(\gamma_T \mid \mathcal{F}_T)\) for the relevant element of \(P_{T\mid T}\), so that \(\sqrt{P_{\gamma,T\mid T}}\) is the conditional standard deviation of the latent slope \(\gamma_T\) given the information available at \(T\). Define the standardised signal \[z_{\gamma,T} = \frac{\hat\gamma_{T\mid T}}{\sqrt{P_{\gamma,T\mid T}}}.\] The reinitialization warning is triggered when \[z_{\gamma,T} > c_{reinit}, \qquad c_{reinit} = 2,\] that is, when \(\hat\gamma_{T\mid T}\) exceeds zero by approximately two conditional standard deviations. This threshold is an operational rule of thumb rather than a fixed statistical decision rule, and can be adjusted depending on the noisiness of the data and the desired trade-off between early warning and false reinitialization.

Let \(d\) denote the first forecast origin at which this trigger condition holds. The reset date \(r\) is then defined as the start of the current positive spell in the smoothed slope estimate computed on the sample ending at \(d\): \[r = \min\{t \le d : \hat\gamma^S_{u\mid d} > 0 \text{ for all } u=t,\ldots,d\},\] where \(\hat\gamma^S_{u\mid d}\) denotes the smoothed estimate of the slope at date \(u\), using observations only up to and including date \(d\). Thus the trigger date \(d\) and the reset date \(r\) are distinct: \(d\) is the date on which the warning first fires in real time, whereas \(r\) is the earlier date at which the cumulative series is reinitialized. In the notation below, \(t=r\) denotes the reset date and \(r_0\) the date at which the cumulative series is set to 0. Then:

\[\begin{align} \ln g_t &= \ln y_t - \ln Y_{t-1} & t=1, \ldots, r & \notag \\ \ln g_t^r &= \ln y_t - \ln Y_{t-1}^r & t=r+1, \ldots, T & (10)\\ Y_{t}^{r}&=Y_{t-1}^{r}+y_{t} & t=r,\ldots,T & (11) \end{align}\]

where \(Y_{t}^{r}\) denotes the cumulative cases after re-initialization. We set \(Y_{r-1}^{r}=0\), so that the growth rate of cumulative cases is available from \(t=r+1\) onwards. Note that \(Y_{t}^r = Y_{t} - Y_{r_0}\).

The gap between the two series becomes apparent by writing

\[\begin{equation} \ln g_t^r = \ln g_t + \ln \frac{Y_{t-1}}{Y_{t-1}^r} = \ln g_t + \ln \frac{Y_{t-1}}{Y_{t-1} - Y_{r_0}} \;\;\;\; t=r+1, \ldots, T \;\;\;\;\;\;\;\; (12) \end{equation}\]

In the next section, where we illustrate the working of the program, it can be seen that in contrast to the original \(\ln g_t\) series, which continues to increase, the reinitialized \(\ln g_t\) series begins to decrease from the reinitialization date. The reinitialization enforces the canonical Gompertz curve with the log of growth rate of cumulative cases sloping down.

Reinitializing the model

We reinitialize the model by specifying the appropriate prior distribution for the initial states: \(\alpha^r_1 \sim N(a^r_1,P^r_1)\) with \(a_1^r = (a_{\delta,1}^r,a_{\gamma,1}^r)'\) and \(P_1^r\) defined as follows. Let \(\mathcal{F}_t = \ln g_t, \ln g_{t-1}, \ldots, \ln g_1\) and define \(a_t = E(\alpha_t|\mathcal{F}_{t-1})\), \(a_{t|t} = E(\alpha_t|\mathcal{F}_t)\), \(P_t = Var(\alpha_t|\mathcal{F}_{t-1})\), \(P_{t|t} = Var(\alpha_t|\mathcal{F}_t).\) Then,

\[\begin{align*} a^r_{\delta,1} &= a_{\delta,r+1} + \ln(Y_r/y_r) \\ a^r_{\gamma,1} &= 0\\ P^r_1 &= P_{r+1}, \end{align*}\]

where \(a_{\delta,r+1}\) and \(P_{r+1}\) are obtained from the non-reinitialized model estimated over \(t=1, \ldots, r\) via the usual Kalman filter recursions. Adding \(\ln(Y_r/y_r)\) to \(a_{\delta,r+1}\) corrects for the shift down in the log cumulative cases caused by reinitializing the cumulative case series. Setting \(a^r_{\gamma,1} = 0\) ensures the model starts off with exponential, rather than super-exponential, growth.

We reinitialize the model through the priors in this way rather than simply re-estimating the model from scratch for two reasons. First, it allows us to impose a proper (rather than diffuse) prior centered on zero for \(\gamma\), so that the starting point is exponential growth. Second, it enables us to make use of data from before the reinitialization date. One needs a reasonable sample size for the estimated model and forecasts to be reliable, but if a new wave is taking off, forecasts need to be generated quickly. This was particularly true with the emergence of the Omicron variant, which caused an explosive increase in infection over a short period of time.

We do not re-estimate the \(\sigma^2_{\varepsilon}\) or \(\sigma^2_{\zeta}\) parameter in the reinitialized model. Rather, we use the values estimated in the original model over \(t=1, \ldots, r\). The one-step-ahead prediction error at \(t=r\) is the same in both the initialized and reinitialized models, but after \(t=r\), the prediction errors diverge.

The reinitialization procedure is very similar in the case where we have seasonal terms. If we let \(\alpha_{s,t}\) be the vector of seasonal states and maintain an analogous notation to that above, the prior mean of the seasonal components in the reinitialized model is \[a^r_{s,1} = a_{s,r+1}.\] The prior variance of \(\alpha^r_1\) remains \(P_{r+1}\) where \(P_{r+1}\) is appropriately re-defined to include the seasonal term, as described in the Appendix.

Leading indicator model

The basic dynamic Gompertz model uses only the cumulative target series. A useful extension is the leading-indicator model, in which a second series is used to improve forecasts of an outcome of interest. The motivation is simple: many indicators move earlier than the target, so their information can be exploited to improve short-horizon forecasts. In the epidemic setting, for example, hospital admissions often respond with a lag to the underlying growth in cases, while in marketing applications early sales of related products may foreshadow demand for a new item. The package therefore allows the target series to be modeled jointly with a lagged auxiliary series rather than treating the target in isolation.

Let \(\ln g_{1,t-k}\) denote the log growth rate of the leading series and \(\ln g_{2,t}\) the log growth rate of the target series. The leading-indicator specification in tsgc is

\[\begin{align} \ln g_{1,t-k} &= \delta_t + \psi_t + \varepsilon_{1t}, & t = k+1, \ldots, T+k, & (16)\\ \ln g_{2,t} &= \delta_t + \varepsilon_{2t}, & t = k+1, \ldots, T, & (17) \end{align}\]

where \(k\) is the user-specified lag, \(\delta_t\) is the common trend shared by both series and evolves exactly as in the dynamic Gompertz model (equations (3)-(4)), \(\psi_t\) is a target-independent deviation specific to the leading equation, and \(\varepsilon_{1t}\) and \(\varepsilon_{2t}\) are mutually independent Gaussian observation disturbances. The lag \(k\) is an alignment parameter: after shifting the leading series by \(k\) periods, the lead and target growth rates are treated as two noisy measurements of a common latent growth component, rather than the leading series entering as an exogenous regressor. Estimation is by maximum likelihood using the Kalman filter, exactly as for the univariate model, and the resulting object is a FilterResultsLI object, which supports the same forecasting, plotting, and diagnostic methods as the FilterResults class used for the dynamic Gompertz model.

The role of n.lag is therefore to align the information in the leading series with the target series, in integer positions rather than calendar time. If the indicator is known to anticipate the target by several periods, a larger lag is appropriate. The package also allows the user to supply exogenous covariates through xpred_lead and xpred_targ at estimation time, and their future values by setting the xpred_lead.new and xpred_targ.new fields directly on the fitted FilterResultsLI object before forecasting, making the model suitable for settings where both the lead and target series are influenced by observed external drivers.

The empirical benefit of this extension is that it can improve short-horizon predictions when the relationship between the two series is stable. In practice, one usually compares several lag values and evaluates them using a rolling-origin experiment, because the appropriate lead-lag relationship is seldom known in advance. Worked examples of the leading-indicator model, including the England cases-hospitalisations application, a comparison against the dynamic Gompertz model, lag selection by rolling cross-validation, and extensions to other data frequencies, are given in the Illustration section below.

Setting up the forecasting exercise

We begin by specifying a number of options for the forecasting exercise, as defined below.

• Y is the data, in the form of an idx_series of cumulative confirmed cases. In this example the object holding this series is called gauteng_idx, with gauteng_cal the accompanying idx_calendar.

• estimation.pos.start is the integer position of the first observation in the sample to be used for estimating the model, obtained here from a calendar date via idx_to_pos(). By default, if start is left unspecified when constructing the model, it is the first position in Y.

• estimation.pos.end is the integer position of the last observation in the sample to be used for estimating the model, again obtained via idx_to_pos(). By default, if end is left unspecified, it is the last position in Y.

• n.forecasts is the number of days or periods for which forecasts are to be made. E.g., if n.forecasts = 14, forecasts will be generated for up to 14 positions following estimation.pos.end.

• q is the signal-to-noise ratio, which controls the smoothness of the estimated trend. A lower value will lead to more smoothness. By default, we use q = NULL, so that it is estimated. Alternatively, we can set q. For example, in our experience, q = 0.005 ensures a good balance between the smoothness of the trend and the speed with which changes in estimates respond to new observations.

• confidence.level sets the coverage used to construct the state-implied forecast bands for \(\ln(g_t)\), which are then used to generate the forecast bands for the forecasts themselves. Here, we use 0.68, corresponding to the probability that the forecast lies within one standard deviation of the point forecast under the (approximating) Normal distribution used to construct the bands – see the note on interpretation in the Theory section above.

• plt.length sets a truncation position to enhance the clarity of plots, e.g. showing only the last 30 positions of the estimation sample. The range for plotting can be set as plt.length positions up to estimation.pos.end.

In this example the data is the cumulative confirmed cases time series for Gauteng. The start and end positions (estimation.pos.start and estimation.pos.end) that define the sample used for estimation are chosen as appropriate for the exercise, translated from calendar dates via idx_to_pos(). In this example, we begin with the sample period set from 1 February to 3 May 2021. This marks the beginning of the third wave in Gauteng as can be seen in Figure 1. The options are specified as below.

Y <- gauteng_idx
cal <- gauteng_cal

estimation.pos.start <- idx_to_pos(cal, "2021-02-01")
estimation.pos.end   <- idx_to_pos(cal, "2021-05-03")
n.forecasts <- 14
q <- 0.005
confidence.level <- 0.68
plt.length <- 30

Estimation

We begin by selecting the data series (Y) for the defined sample period.

The estimation period is specified directly when creating the SSModelDynamicGompertz object using the start and end arguments (integer positions in Y), so no manual subsetting of the data is required. Passing calendar = cal attaches the idx_calendar to the resulting model and, in turn, to its FilterResults output, so that later plots are automatically labelled with real dates rather than raw positions.

We then estimate the model using a diffuse prior distribution for the initial state vector. The signal-to-noise ratio can be left as a free parameter to be estimated, as in the code below.

model_q <- SSModelDynamicGompertz(Y = Y, sea.period = 7,
                                  start = estimation.pos.start,
                                  end = estimation.pos.end,
                                  calendar = cal)
res_q <- estimate(model_q)

In the rest of this example we estimate the model setting the signal-to-noise ratio at \(0.005.\) As mentioned, in our experience this value strikes a useful balance between the smoothness of the estimate of the slope parameter \(\gamma\), and the speed with which it adapts to new observations.

model <- SSModelDynamicGompertz(Y = Y, q = q, sea.period = 7,
                                start = estimation.pos.start,
                                end = estimation.pos.end,
                                calendar = cal)
res <- estimate(model)

Model diagnostics

Before looking at forecasts, it is worth checking that the model actually fitted sensibly. summary() on a FilterResults/ FilterResultsLI object reports the estimation window, model states and their standard errors, and the estimated variance parameters (including the signal-to-noise ratio); tsgc::print_model_diagnostics() supplements this with the fitted log-likelihood and simple recursive-residual diagnostics.

summary(res)

## Summary of FilterResults Object
## Model Details:
##   - Estimation start position: 329
##   - Estimation end position: 420
##   - Model States and Standard Errors
## Smoothed values of states and standard errors at time n = 92:
##             Estimate   Std. Error
## level       -6.854409   0.114903 
## slope        0.059714   0.033423 
## sea_trig1   -0.085890   0.030406 
## sea_trig*1   0.466831   0.030605 
## sea_trig2    0.206475   0.030267 
## sea_trig*2   0.173606   0.030246 
## sea_trig3    0.148638   0.030259 
## sea_trig*3   0.000876   0.030200 
## 
##   - Variance parameter estimates
## Observation equation noise: 0.04129
## State transition equation noise: 0.0002064
## Signal-to-Noise Ratio (q): 0.005
## Seasonality noise: 1.333e-07

tsgc::print_model_diagnostics(res)

## ---- Model diagnostics ----
##   Log-likelihood: -15.6996 
##   Recursive residuals: mean = 0.1073 , sd = 1.179 
## ----------------------------

Results

We can now plot the forecast of \(\ln g_t\) – the log of the growth rate of \(Y\), the cumulative cases – which is the transformation of the series that is taken to the model, and we can compare these forecasts to the actual \(\ln g_t\) series. We do this by passing the output (res) of the estimation step along with an evaluation sample to a plotting function tsgc::plot_log_forecast.

tsgc::plot_log_forecast(
  res = res,
  Y = Y,
  n.ahead = n.forecasts,
  plt.start = tail(res$index, 1) - plt.length
)

Figure 2: Fourteen-day forecast of \(ln(g_t)\) from 4 May 2021, for Gauteng province in South Africa.

From these results we can recover the forecasts of new cases from 4 May 2021, with their forecast bands.

tsgc::plot_forecast(
  res = res,
  n.ahead = n.forecasts,
  confidence.level = confidence.level,
  plt.start = tail(res$index, 1) - plt.length,
  series.name = "cases"
)

Figure 3: Fourteen-day forecast of new cases from 4 May 2021 for Gauteng province in South Africa. The shaded region is a state-implied forecast band obtained by propagating uncertainty in the latent log growth rate.

To assess accuracy, we plot these forecasts against the actual new cases, that have been held back from the estimation sample, using the plot_holdout function. The model forecasts are compared with the first differences of Y, the cumulative series, over the forecast window.

# Update estimation end position
estimation.pos.end <- idx_to_pos(cal, "2021-04-19")

# Reestimate the model
model <- SSModelDynamicGompertz$new(Y = Y, q = q, sea.period = 7,
                                    start = estimation.pos.start,
                                    end = estimation.pos.end,
                                    calendar = cal)
res_eval <- estimate(model)

tsgc::plot_holdout(
  res = res_eval,
  Y = Y,
  n.ahead = 14,
  confidence.level = 0.68,
  series.name = "cases"
)

Figure 4: Accuracy of the fourteen-day forecast of new cases from 20 April 2021 for Gauteng province in South Africa.

Figure 4 shows that the forecasts were accurate over the first seven days, with a mean absolute percentage error (MAPE) of 13.9%. Note that reported cases were unusually low on 27 April due to the fact that it is Freedom Day and a public holiday in South Africa. That day aside, over the six days from 28 April the MAPE was 12.8%. Over the full 14 days of the forecasts, the MAPE was 27%.

As discussed earlier, the reproduction numbers \(R_t\) and their 68% approximate pointwise uncertainty bands can be calculated using estimate_r0(). It returns a data frame containing the integer position, fitted \(R_t\), and lower and upper bounds. When the fitted result has a calendar, as here, the corresponding date is included as well. The estimates can be plotted directly with plot_r0(), as shown in Figure 5.

gen_int <- 4
r.t <- estimate_r0(res, gen_int, n.ahead = 7)
r.t

##         Date Position       fit     lower     upper
## 1 2021-04-27      414 0.8309205 0.7274680 0.9490848
## 2 2021-04-28      415 0.9449307 0.8272860 1.0793050
## 3 2021-04-29      416 1.0433477 0.9134504 1.1917171
## 4 2021-04-30      417 1.1241605 0.9842048 1.2840181
## 5 2021-05-01      418 1.1274117 0.9870550 1.2877267
## 6 2021-05-02      419 1.1937255 1.0451162 1.3634661
## 7 2021-05-03      420 1.2751620 1.1164161 1.4564804

plot_r0(res, gen_int = gen_int, n.ahead = 7,
        title = "Reproduction numbers for the 7-day period to 3 May 2021")

Figure 5: Reproduction numbers for the 7-day period to 3 May 2021, for Gauteng province in South Africa.

The dated data frame printed above also makes the \(R_t\) values and their confidence interval available directly, without having to read them from the plot.

A CSV file (y-forecast) is written to the directory specified. The forecast options specified earlier are retained.

res.dir <- here::here(here::here(), "results")

tsgc::write_results(
 res = res,
 res.dir = res.dir,
 n.ahead = n.forecasts,
 confidence.level = confidence.level
)

Reinitialization

In all countries, the coronavirus pandemic was characterized by a series of recurring waves due to a combination of biological, behavioral, and environmental reasons. In an epidemic, the onset of a new wave is signalled when the slope parameter \(\gamma\), which measures of the growth rate of the growth rate of new cases, rises above zero for a sustained period. Such a super-exponential phase of the epidemic in which the growth rate of new cases is itself increasing over time is typically short.

This section illustrates the reinitialization procedure which allows us to apply the model to the new wave as it begins, without jettisoning information from the wave that has just ended. We extend the estimation window to 25 June 2021, by which date the third wave is well on course with its peak within sight (see Figure 1). All other options remain the same.

estimation.pos.end <- idx_to_pos(cal, "2021-06-25")

Triggering reinitialization

As set out above, reinitialization is triggered at the first forecast origin \(T\) at which the standardised signal \(z_{\gamma,T} = \hat\gamma_{T\mid T}/\sqrt{P_{\gamma,T\mid T}}\) exceeds the threshold \(c_{reinit}=2\), that is, when the estimate of the slope parameter, \(\gamma_t\), breaches two standard errors above zero. Reinitialization is then backdated to \(r\), the start of the current positive spell in the smoothed slope estimate computed on the sample ending at the trigger position. In applying this rule there is a choice between the filtered slope estimate and the smoothed slope estimate. Experimental results suggest that the greater noise in the filtered estimate of \(\gamma_t\) often triggers reinitialization too early. The smoothed estimate is more reliable, and is what we use to compute \(z_{\gamma,T}\) below.

Figure 6 shows that for the third wave in Gauteng, the smoothed slope estimate exceeded twice its standard error on 1 May 2021, having risen above zero on 21 April 2021, giving trigger date \(d =\) 1 May 2021 and reset date \(r =\) 21 April 2021.

# Re-estimate model
model <- SSModelDynamicGompertz$new(Y = Y, q = q, sea.period = 7,
                                    start = estimation.pos.start,
                                    end = estimation.pos.end,
                                    calendar = cal)
res <- estimate(model)

# Extract the smoothed slope and its standard deviation as idx_series,
# then convert to a data.frame with real dates for plotting.
smoothed.slope.full <- idx_series(res$output$alphahat[, "slope"], start = res$index[1])
V.smoothed <- get_V(res$output)
i.slope <- grep("slope", colnames(res$output$alphahat))
smoothed.P.slope <- idx_series(V.smoothed[i.slope, i.slope, ], start = res$index[1])

common_pos <- idx_positions(smoothed.slope.full)
d2 <- idx_series(
  cbind(
    smoothed.slope = idx_values(smoothed.slope.full),
    sd.smoothed.slope = sqrt(idx_values(smoothed.P.slope)),
    sd.smoothed.slope.1.5 = 1.5 * sqrt(idx_values(smoothed.P.slope)),
    sd.smoothed.slope.2 = 2 * sqrt(idx_values(smoothed.P.slope))
  ),
  start = common_pos[1]
)

d2.mat <- as.matrix(idx_values(d2))
d2.df <- data.frame(
  Date = idx_to_date(cal, idx_positions(d2)),
  smoothed.slope = d2.mat[, "smoothed.slope"],
  sd.smoothed.slope = d2.mat[, "sd.smoothed.slope"],
  sd.smoothed.slope.1.5 = d2.mat[, "sd.smoothed.slope.1.5"],
  sd.smoothed.slope.2 = d2.mat[, "sd.smoothed.slope.2"]
)

d2.df <- d2.df[d2.df$Date >= as.Date("2021-02-15"), ]

# zt = smoothed slope minus its own two-standard-error threshold (the
# lower bound of the two-SE interval).
d2.df$zt <- d2.df$smoothed.slope - d2.df$sd.smoothed.slope.2

trigger.df <- d2.df %>%
  mutate(prev_zt = lag(zt)) %>%
  filter(zt > 0 & prev_zt <= 0)
# Triggered when the two-SE lower bound of the smoothed slope crosses zero

reinit_zero.df <- d2.df %>%
  mutate(prev_smoothed.slope = dplyr::lag(smoothed.slope)) %>%
  filter(Date < min(trigger.df$Date) &
           (smoothed.slope > 0 & prev_smoothed.slope < 0)) %>%
  arrange(desc(Date)) %>%
  slice(1)
# Reinitialisation on April 21.

# Create plot
ggplot(data = d2.df,
      aes(x = Date)) +
  geom_line(aes(y = smoothed.slope,
                color = "smoothed.slope"), linewidth=.5) +
  geom_line(aes(y = sd.smoothed.slope,
                color = "sd.smoothed.slope"),
            linetype = "solid",
            linewidth=.25) +
  geom_line(aes(y = sd.smoothed.slope.1.5,
                color = "sd.smoothed.slope.1.5"),
            linetype = "solid",
            linewidth=.25) +
  geom_line(aes(y = sd.smoothed.slope.2,
                color = "sd.smoothed.slope.2"),
            linetype = "solid",
            linewidth=.5) +
  scale_y_continuous(n.breaks = 10) +
  geom_hline(yintercept = 0,
             linetype = "solid",
             color = "black",linewidth = 1)+
  geom_vline(data = trigger.df,
             aes(xintercept = Date),
             linetype = "solid", linewidth = .5,
             color = "black")+
  ylab("Slope") +
  scale_x_date(date_breaks = "10 days") +
  scale_color_manual(name='',
                     values=c('smoothed.slope'='red',
                              'sd.smoothed.slope'='blue',
                              'sd.smoothed.slope.1.5'='green',
                              'sd.smoothed.slope.2'='black'))+
  theme_light(base_size = 11)+
  #theme(legend.position = c(.88, 0.25))+
  theme(legend.title=element_text(size=2),
        legend.text=element_text(size=6))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
        plot.title = element_text(face = "bold"))

Figure 6: Trigger for reinitialization: reinitialization is triggered when the two-standard-error lower threshold of the smoothed slope (slope - 2*se) crosses zero. Lines shown are k-SE thresholds, not confidence bands around the slope. Reset date selected from a single retrospective episode; not a validated general rule.

The position for reinitialization is set accordingly, translating the reset date to an integer position via idx_to_pos().

reinit.pos <- idx_to_pos(cal, "2021-04-21")

Estimating the reinitialized model

A reinitialized SSModelDynamicGompertz uses the same arguments as the non-reinitialized model, with the addition of the reinit.idx argument, an integer position.

model <- SSModelDynamicGompertz$new(
  Y = Y, q = q, sea.period = 7,
  start = estimation.pos.start,
  end = estimation.pos.end,
  reinit.idx = reinit.pos,
  calendar = cal
)
res.reinit <- estimate(model)

We generate the reinitialized data series by setting cumulative cases to \(0\) at the appropriate point, as discussed in the Theory section, and extract the evaluation sample from the reinitialized series as below.

Estimating the model with the reinitialized series, the actual and forecast \(\ln(g_t)\) can be plotted as in Figure 7, which is analogous to Figure 2 for the non-reinitialized series.

tsgc::plot_log_forecast(
    res = res.reinit,
    Y = Y,
    n.ahead = n.forecasts,
    plt.start = tail(res.reinit$index, 1) - plt.length,
    title='Forecast of log growth rate after reinitialization.'
)

Figure 7: Forecast of \(ln(g_t)\) after reinitialization.

Comparing forecasts, without reinitialization the 14-day forecast MAPE is 41.9% (see Figure 8). With reinitialization, it falls to 20.2% (see Figure 9). The corresponding MAPE figures for 7-day forecasts is 15.2% without reinitialization and 9.5% with reinitialization.

tsgc::plot_holdout(
  res = res,
  Y = Y,
  n.ahead = 14,
  confidence.level = 0.68,
  series.name = "cases"
)

Figure 8: Forecast accuracy of the model without reinitialization over the hold-out sample period: 14 days from 25 June 2021.

tsgc::plot_holdout(
  res = res.reinit,
  Y = Y,
  n.ahead = 14,
  confidence.level = 0.68,
  series.name = "cases"
)

Figure 9: Forecast accuracy of the model with reinitialization over the hold-out sample period: 14 days from 25 June 2021.

We can also compare the standard (no reinitialization) and reinitalized forecasts in a single plot:

tsgc::plot_compare_forecast(
  results = list(res.reinit, res),
  n.ahead = 14,
  actual = Y
)

Figure 10: Comparison of standard (green) and reinitialized (blue) forecasts

Just as for the standard (non-reinitialized) model, the returned estimation results contained in res.reinit are a FilterResults object, and can be written to CSV using the write_results function.

Leading indicator model: England cases and hospitalisations

We now illustrate the leading-indicator model described above. We load the England data, which contains daily cumulative COVID-19 cases and hospital admissions, convert it to idx_series/idx_calendar form, and specify an estimation window running from 30 April 2021 to 24 July 2021, with cases as a four-day leading indicator for hospitalisations.

data(england, package = "tsgc")

conv <- xts_to_idx(england[, 1:2])
eng <- conv$series
eng_cal <- conv$calendar

est.start.eng <- idx_to_pos(eng_cal, "2021-04-30")
est.end.eng   <- idx_to_pos(eng_cal, "2021-07-24")
n.lag.eng     <- 4
n.forecasts.eng <- 7
plt.length.eng  <- 14

A preliminary plot of the two series on the log scale is useful for inspecting comovement before formal estimation.

mod.eng.plot <- SSModelLeadingIndicator$new(eng, n.lag = n.lag.eng, calendar = eng_cal)
plot(
  mod.eng.plot,
  title = "Daily COVID cases and hospitalisations (England)",
  series.name.lead = "Cases",
  series.name.target = "Hospitalisations",
  take.log = TRUE
)

Figure 11: Daily COVID cases and hospitalisations, England.

We then estimate the baseline leading-indicator model, with cases as the lead series (LeadIndCol = 1) and weekly seasonality.

mod.eng <- SSModelLeadingIndicator$new(
  Y = eng,
  n.lag = n.lag.eng,
  LeadIndCol = 1,
  sea.period = 7,
  start = est.start.eng,
  end = est.end.eng,
  calendar = eng_cal
)
res.eng <- estimate(mod.eng)

FilterResultsLI’s summary() prints the same kind of diagnostics as FilterResults$summary() above, but its observation-equation noise H is a \(2\times2\) matrix (one variance per equation, for the lead and target series) rather than a single scalar; it is printed as a matrix accordingly.

summary(res.eng)

## Summary of FilterResultsLI Object
## Model Details:
##   - Estimation start position: 407
##   - Estimation end position: 492
##   - Model States and Standard Errors
## Smoothed values of states and standard errors at time n = 86:
##                     Estimate   Std. Error
## level               -6.227326   0.039235 
## slope                0.072775   0.184027 
## level.LDLlead        1.323149   0.068104 
## sea_trig1.LDLlead    0.161333   0.041468 
## sea_trig*1.LDLlead   0.103834   0.041468 
## sea_trig2.LDLlead   -0.019902   0.019065 
## sea_trig*2.LDLlead  -0.074915   0.019092 
## sea_trig3.LDLlead   -0.039680   0.016282 
## sea_trig*3.LDLlead   0.008766   0.016482 
## sea_trig1.LDLtarg   -0.093289   0.036842 
## sea_trig*1.LDLtarg   0.004641   0.036842 
## sea_trig2.LDLtarg   -0.019120   0.011347 
## sea_trig*2.LDLtarg   0.009455   0.011347 
## sea_trig3.LDLtarg   -0.006348   0.007297 
## sea_trig*3.LDLtarg  -0.010262   0.007297 
## 
##   - Variance parameter estimates
## Observation equation noise:
##        [,1] [,2]
## [1,] 0.0077    0
## [2,] 0.0000    0
## State transition equation noise: 0.03233
## Seasonality noise: 0.005768

tsgc::print_model_diagnostics(res.eng)

## ---- Model diagnostics ----
##   Log-likelihood: 42.3087 
##   Recursive residuals: mean = -0.03584 , sd = 1.236 
## ----------------------------

The fitted object res.eng can be used to generate the same forecast and holdout panels as for the dynamic Gompertz model.

tsgc::plot_log_forecast(
  res.eng,
  Y = eng,
  n.ahead = n.forecasts.eng,
  plt.start = est.end.eng - plt.length.eng,
  title = "Forecast of log growth rate of hospitalisations (England)"
)

Figure 12: Forecast of log growth rate of hospitalisations, England.

tsgc::plot_forecast(
  res.eng,
  n.ahead = n.forecasts.eng,
  plt.start = est.end.eng - plt.length.eng,
  series.name = "Hospitalisations",
  title = "Forecast of hospitalisations (England)"
)

Figure 13: Forecast of hospitalisations, England.

tsgc::plot_holdout(
  res.eng,
  Y = eng,
  n.ahead = n.forecasts.eng,
  series.name = "Hospitalisations",
  title = "Accuracy: forecast of hospitalisations (England)"
)

Figure 14: Accuracy of the seven-day forecast of hospitalisations, England.

Adding regressors and future covariate paths

The same specification can be enriched with exogenous predictors. Here we add contemporaneous weather variables (wind speed, relative humidity, temperature, and maximum temperature) to both the lead and target equations, and then supply their future values before forecasting.

data(england_weather_2021, package = "tsgc")

conv <- xts_to_idx(england_weather_2021[, 1:4],
                    start.pos = idx_to_pos(eng_cal, zoo::index(england_weather_2021)[1]))
england_weather_idx <- conv$series

xpred_lead <- xpred_targ <- england_weather_idx

mod.eng.x <- SSModelLeadingIndicator$new(
  eng,
  n.lag = n.lag.eng,
  xpred_lead = xpred_lead,
  xpred_targ = xpred_targ,
  sea.period = 7,
  start = est.start.eng,
  end = est.end.eng,
  calendar = eng_cal
)
res.eng.x <- estimate(mod.eng.x)
res.eng.x$xpred_lead.new <- england_weather_idx
res.eng.x$xpred_targ.new <- england_weather_idx

tsgc::plot_forecast(
  res.eng.x,
  n.ahead = n.forecasts.eng,
  plt.start = est.end.eng - plt.length.eng,
  title = "Forecast of hospitalisations\nwith regressors (weather, oracle/realised)",
  series.name = "Hospitalisations"
)

Figure 15: Forecast of hospitalisations with weather regressors (weather, oracle/realised), England.

In practice, the recommended sequence is to establish a sensible lag and seasonal specification first, inspect the resulting forecasts, and only then augment the model with regressors if additional covariates are available and future paths can be supplied.

Comparing the dynamic Gompertz and leading-indicator models

We now compare a UK-only dynamic Gompertz model with a leading-indicator model that uses Italian case growth as a lead for UK case growth, using daily COVID-19 cases in Italy and the UK from the first-wave period. Italy experienced the first large European outbreak earlier than the UK, making it a natural candidate for a short-horizon lead series during the first wave.

data(ukitaly, package = "tsgc")

conv <- xts_to_idx(ukitaly)
ukitaly_idx <- conv$series
ukitaly_cal <- conv$calendar

n.forecasts.uk <- 14
plt.length.uk  <- 30
q.uk   <- 0.005

est.start.uk <- idx_to_pos(ukitaly_cal, "2020-02-25")
est.end.uk   <- idx_to_pos(ukitaly_cal, "2020-04-01")
Yuk          <- idx_series(idx_values(ukitaly_idx)[, "UK"], start = ukitaly_idx$start)

mod.ukit.plot <- SSModelLeadingIndicator$new(ukitaly_idx, n.lag = 4, calendar = ukitaly_cal)
plot(
  mod.ukit.plot,
  title = "Daily COVID cases in UK and Italy",
  series.name.lead = "Italy",
  series.name.target = "UK",
  take.log = FALSE
)

Figure 16: Daily COVID cases in UK and Italy.

The UK-only dynamic Gompertz model is estimated on UK daily cases alone, and the leading-indicator model uses Italy as the lead series with a 14-day lag, motivated by the observed timing of the first-wave evolution in Italy relative to the UK. Both models are estimated with weekly seasonality (sea.period = 7, the package default for daily data), specified explicitly here so that the two specifications are directly comparable; a fuller factorial comparison across seasonal and other tuning choices is beyond the scope of this vignette.

mod.uk.gomp <- SSModelDynamicGompertz$new(
  Y = Yuk, q = q.uk, sea.period = 7,
  start = est.start.uk, end = est.end.uk,
  calendar = ukitaly_cal
)
res.uk.gomp <- estimate(mod.uk.gomp)

n.lag.uk <- 14
mod.uk.lead <- SSModelLeadingIndicator$new(
  Y = ukitaly_idx, n.lag = n.lag.uk, sea.period = 7,
  start = est.start.uk, end = est.end.uk,
  calendar = ukitaly_cal
)
res.uk.lead <- estimate(mod.uk.lead)

tsgc::plot_compare_forecast(
  list(res.uk.gomp, res.uk.lead),
  actual = Yuk,
  n.ahead = n.forecasts.uk
)

Figure 17: Comparison of dynamic Gompertz and leading-indicator forecasts, UK daily cases.

tsgc::plot_holdout(
  res.uk.lead, Y = ukitaly_idx, n.ahead = n.forecasts.uk,
  title = "Accuracy: leading-indicator forecast of UK cases",
  series.name = "UK cases"
)

Figure 18: Accuracy of the leading-indicator forecast, UK daily cases.

In this first-wave window, the leading-indicator model gives substantially smaller holdout errors than the UK-only dynamic Gompertz model.

est.end.cv.uk <- idx_to_pos(ukitaly_cal, "2020-04-15")

cv_models <- list()

cv_models[["Naive_last_value"]] <- SSModelDynamicGompertz$new(
  Y = Yuk, q = 0, sea.period = 7,
  start = est.start.uk, end = est.end.cv.uk,
  calendar = ukitaly_cal
)

cv_models[["RW_growth"]] <- SSModelDynamicGompertz$new(
  Y = Yuk, q = q.uk, sea.period = 0,
  start = est.start.uk, end = est.end.cv.uk,
  calendar = ukitaly_cal
)

cv_models[["Vanilla_q"]] <- SSModelDynamicGompertz$new(
  Y = Yuk, q = q.uk, sea.period = 7,
  start = est.start.uk, end = est.end.cv.uk,
  calendar = ukitaly_cal
)

cv_models[["Vanilla_ar1"]] <- SSModelDynamicGompertz$new(
  Y = Yuk, sea.period = 7,
  start = est.start.uk, end = est.end.cv.uk, ar1 = TRUE,
  calendar = ukitaly_cal
)

for (i in c(7, 10, 14, 16, 18, 21)) {
  cv_models[[paste0("Lag", i)]] <- SSModelLeadingIndicator$new(
    Y = ukitaly_idx, sea.period = 7,
    start = est.start.uk, end = est.end.cv.uk, n.lag = i,
    calendar = ukitaly_cal
  )
}

n.ahead.cv  <- 5
gap.cv      <- n.ahead.cv   # non-overlapping horizons
n.select.cv <- 2
n.report.cv <- 1

report.end.cv <- est.end.cv.uk - (n.report.cv - 1) * gap.cv - n.ahead.cv
# The last SELECTION origin is select.end.cv + (n.select.cv-1)*gap.cv,
# with its forecast horizon reaching + n.ahead.cv further. Requiring
# that to finish strictly before the REPORTING origin (report.end.cv)
# means subtracting n.select.cv * gap.cv.
select.end.cv <- report.end.cv - n.select.cv * gap.cv - n.ahead.cv
stopifnot(select.end.cv > est.start.uk)

# Explicit disjointness check: the last SELECTION fold's forecast
# horizon must end strictly before the REPORTING block's origin.
last_selection_origin      <- select.end.cv + (n.select.cv - 1) * gap.cv
last_selection_horizon_end <- last_selection_origin + n.ahead.cv
stopifnot(last_selection_horizon_end < report.end.cv)

# Feasibility check for the highest lag tested: SSModelLeadingIndicator's
# estimate() needs at least n.lag + 2 usable positions from est.start.uk
# before the SELECTION block's origin (n.lag positions consumed by the
# forward lag-shift, plus at least 1-2 more lost to differencing/NA
# trimming).
max_lag_cv <- 21
stopifnot(select.end.cv - est.start.uk > max_lag_cv + 1)

cv_selection <- cross_val(
  Y = ukitaly_idx,
  model_list = cv_models,
  est.end = select.end.cv,
  criterion = "smape",
  n.ahead = n.ahead.cv,
  n.estimate = n.select.cv,
  gap = gap.cv
)
cv_selection

##               Model    31    36
## 1  Naive_last_value 11.37 20.21
## 2         RW_growth  8.45  8.17
## 3         Vanilla_q  5.01 13.41
## 4       Vanilla_ar1  3.24 10.86
## 5              Lag7  9.79  6.24
## 6             Lag10  4.79  9.01
## 7             Lag14 31.58  4.51
## 8             Lag16  3.90  4.30
## 9             Lag18  3.59  7.71
## 10            Lag21  6.15  4.05

lag_rows <- cv_selection[grepl("^Lag", cv_selection$Model), ]
lag_rows$mean_smape <- rowMeans(lag_rows[, -1, drop = FALSE])
selected_lag_name <- lag_rows$Model[which.min(lag_rows$mean_smape)]
selected_lag_name

## [1] "Lag16"

cv_report_models <- cv_models[c(
  "Naive_last_value", "RW_growth", "Vanilla_q", "Vanilla_ar1", selected_lag_name
)]
cv_reporting <- cross_val(
  Y = ukitaly_idx,
  model_list = cv_report_models,
  est.end = report.end.cv,
  criterion = "smape",
  n.ahead = n.ahead.cv,
  n.estimate = n.report.cv,
  gap = gap.cv
)
cv_reporting

##              Model    46
## 1 Naive_last_value 27.88
## 2        RW_growth 10.80
## 3        Vanilla_q  4.99
## 4      Vanilla_ar1  7.18
## 5            Lag16  2.16

Extensions: non-epidemiological applications and alternative data frequencies

The package is not restricted to daily epidemic data. The same state-space machinery can be applied to quarterly, monthly, or annual data when the underlying series is cumulative and exhibits saturation; what changes with the sampling frequency is only the interpretation of the integer position, forecast horizon, lag length, and seasonal period, not the underlying estimation procedure. For example, the quarterly Nintendo sales data can be used with a seasonal period of 4 and a lagged leading indicator, while monthly trading-app downloads can be treated in the same way with a seasonal period of 12. In both cases, idx_calendar carries the quarterly/monthly step size and unit, and idx_to_pos() still translates calendar dates to the integer positions that start/end/n.lag require.

data(nintendo_sales, package = "tsgc")

y_q_xts <- nintendo_sales[, c("wii", "switch_all")]

nintendo_idx <- idx_series(zoo::coredata(y_q_xts), start = 1L)
nintendo_cal <- idx_calendar(
  anchor = as.Date(zoo::index(y_q_xts)[1]),
  anchor_pos = 1L,
  amount = 1, unit = "quarters",
  posixct = TRUE
)

n.lag.q <- idx_to_pos(nintendo_cal, "2017-01-01") - idx_to_pos(nintendo_cal, "2006-10-01")

mod.q <- SSModelLeadingIndicator$new(
  Y = nintendo_idx,
  sea.period = 4,
  n.lag = n.lag.q,
  start = idx_to_pos(nintendo_cal, "2017-01-01"),
  end = idx_to_pos(nintendo_cal, "2019-10-01"),
  calendar = nintendo_cal
)
res.q <- estimate(mod.q)

tsgc::plot_log_forecast(
  res.q, Y = nintendo_idx, n.ahead = 8,
  title = "Log forecasts of Switch sales"
)

Figure 19: Log forecasts of quarterly Switch sales, using Wii sales as a leading indicator.

data(etrading_apps, package = "tsgc")

y_m_xts <- etrading_apps[, c("DEGIRO", "AvaTrade")]

etrading_idx <- idx_series(zoo::coredata(y_m_xts), start = 1L)
etrading_cal <- idx_calendar(
  anchor = as.Date(zoo::index(y_m_xts)[1]),
  anchor_pos = 1L,
  amount = 1, unit = "months",
  posixct = TRUE
)

n.lag.m <- idx_to_pos(etrading_cal, "2017-07-01") - idx_to_pos(etrading_cal, "2017-01-01")

mod.m <- SSModelLeadingIndicator$new(
  Y = etrading_idx,
  sea.period = 12,
  n.lag = n.lag.m,
  start = idx_to_pos(etrading_cal, "2017-07-01"),
  end = idx_to_pos(etrading_cal, "2021-02-01"),
  calendar = etrading_cal
)
res.m <- estimate(mod.m)

tsgc::plot_forecast(
  res.m, n.ahead = 4,
  title = "AvaTrade monthly downloads",
  series.name = "Downloads"
)

Figure 20: Forecasts of monthly AvaTrade downloads, using DEGIRO downloads as a leading indicator.

Both examples use the same workflow as the daily illustrations above: select an estimation window, fit the model, choose a forecast horizon, and evaluate the resulting forecasts. No change to the state-space specification or estimation procedure is required when moving between frequencies.

mode = "date" (the default here, since gauteng_cal$posixct = TRUE)

Real calendar dates on the x-axis - this is what axis = NULL already produced in Figure 3.

showcase_axis_plot(idx_axis_opts(mode = "date"))

Real calendar dates (mode = “date”).

mode = "position"

The raw integer idx_series position, ignoring the calendar entirely.

showcase_axis_plot(idx_axis_opts(mode = "position"))

Raw integer positions (mode = “position”).

mode = "steps"

Steps from the calendar’s anchor position (position - anchor_pos). Since gauteng_cal has an anchor_name set ("first recorded case", added earlier when the calendar was built), the axis label uses it automatically.

showcase_axis_plot(idx_axis_opts(mode = "steps"))

Steps from the anchor (mode = “steps”).

mode = "time_since"

The pattern-weighted calendar offset from the anchor, expressed in calendar$units rather than a raw step count - here “days since first recorded case”. For a daily series with no gaps this is numerically identical to "steps", but it diverges for series with an irregular pattern (e.g. business days) or a non-day unit (e.g. quarters, months), and it also works for posixct = TRUE calendars where it can be more readable than a real date.

showcase_axis_plot(idx_axis_opts(mode = "time_since"))

Time since the anchor, in calendar units (mode = “time_since”).

Adding an info box

Setting info_box = TRUE (on any of the modes above) appends a caption summarising the calendar: the anchor (and its name, if set), the step size/unit, and the pattern. pattern_n truncates a long pattern to its first n values; by default the full pattern is shown.

showcase_axis_plot(idx_axis_opts(mode = "steps", info_box = TRUE))

Steps from anchor, with an info box summarising the calendar.

If the info box’s pattern were long (e.g. a business-day calendar with weekends skipped), pattern_n keeps the caption readable:

showcase_axis_plot(idx_axis_opts(mode = "steps", info_box = TRUE, pattern_n = 3))

Info box with pattern truncated to its first 3 values via pattern_n.