--- title: "Forecasting epidemic trajectories: Time Series Growth Curves package `tsgc`" author: "Michael Ashby, Andrew Harvey, Paul Kattuman, Edwin Tang, Guglielmo Secchi, Craig Thamotheram" date: "`r Sys.Date()`" output: rmarkdown::html_vignette bibliography: book.bib vignette: > %\VignetteEncoding{UTF-8} %\VignetteIndexEntry{Forecasting epidemic trajectories: Time Series Growth Curves package `tsgc`} %\VignetteEngine{knitr::rmarkdown} abstract: "This paper documents the Time Series Growth Curves (`tsgc`) package for `R`, which is designed for forecasting epidemic trajectories, including the detection of new waves and turning points. The package implements dynamic Gompertz state-space methods, together with a leading-indicator extension that uses a related series moving ahead of the variable of interest to improve short-horizon forecasts, with lag selection supported by rolling-origin cross-validation. The state-space framework applies equally to daily, weekly, monthly, quarterly, or annual data. Estimation and forecasting are defined purely in terms of integer-position (`idx_series`) data; calendar time is reintroduced only as an optional, purely cosmetic layer for plotting. The package is illustrated with COVID-19 confirmed cases, hospitalisation data, a cross-country leading-indicator comparison, and non-epidemic applications such as the diffusion of new products and the adoption of online services. The models are estimated with Kalman-filter techniques and are suitable for short-horizon forecasting of variables whose cumulative form is subject to an unknown saturation level." editor_options: markdown: wrap: 72 --- ```{r setup, include = FALSE} knitr::opts_chunk$set(message = FALSE, warning = TRUE, fig.width = 7, fig.height = 7) library(tsgc) library(dplyr) library(ggplot2) library(ggthemes) library(zoo) library(xts) theme_set(theme_economist_white(gray_bg = FALSE, base_size = 16)) ``` # Introduction Outbreaks of infectious diseases with epidemic potential require real-time responses by public health authorities. Accurate real-time forecasting of the trajectory of the epidemic over the near future is of great value in this regard. The `R` package, `tsgc`, is intended for use in monitoring and forecasting the progress of an epidemic, including the detection of *new waves* and *turning points*. It develops and implements time series growth curve methods first reported in @Harvey2020 (hereinafter referred to as HK). HK develop a class of time series models for predicting future values of a variable which, when cumulated, is subject to an unknown saturation level. In a single wave of an epidemic, as more and more people get infected, the pool of susceptible individuals dwindles. This results in the decline of new infections, and the cumulative number of infections approaches its saturation level. The model can take account of deviations relative to this canonical trajectory due to changes in social behavior and policy. Models in this family are relevant for many other disciplines, such as marketing (when estimating the demand for new products). While attention here is focused on the spread of epidemics and the applications used for illustration relate to coronavirus, this package is designed with a view to wider applicability. Given the number of different modeling approaches for epidemics, there are many notable packages that can be used for monitoring epidemics. For the most part, these seek to model explicitly the mechanism by which the disease spreads through the population. For example, `EpiModel` [@Jenness2018], `EpiEstim` [@cori_etal], `epinowcast` [@sam_abbott_pietro_monticone_2021] to name a few, can be categorized as belonging to the class of *mechanistic* models in the language of philosophy of science, in that they require structural knowledge of the disease spread mechanism in order to obtain predictions. In contrast, the empirical approach implemented in `tsgc` falls into the class of models described as *phenomenological*. Although it is motivated by the archetypal pattern in the dynamics of disease spread, it does not rely on structural assumptions derived from epidemiological theory. There are advantages to not requiring assumptions about values of parameters relating to, inter alia, disease infectiousness, disease severity, or contact patterns, which are difficult to pin down with sufficient precision, especially in real-time during an epidemic. Our approach makes minimal assumptions and merely requires past observations of the epidemic variable of interest, to which we apply time-series methods to provide predictions over short future time horizons. The model can be estimated quickly and straightforwardly, and subjected to standard diagnostic tests. A statistical model of this type is a useful complement to mechanistic models that attempt to describe the epidemic in terms of underlying processes. Section 2 sets out the state space formulation of the dynamic Gompertz growth curve and the way nowcasts and forecasts are obtained from predictive recursions. It is then shown how these numbers translate into estimates of the instantaneous reproduction number $R_t$. Section 3 explains how multiple waves can be accommodated by reinitializing the series at the start of new waves. The start of a new wave is not obvious in real time but a rule for triggering reinitialization that works well in practice is presented. Section 4 describes the functionality of `tsgc`, including the leading-indicator extension. Section 5 sets out a full working example of the package, beginning with COVID-19 data for Gauteng province in South Africa and continuing with the leading-indicator model, a comparison against the dynamic Gompertz model, lag selection by rolling cross-validation, the forecast accuracy measures used to evaluate the package's models, and the package's applicability to data of non-daily frequency and non-epidemological nature. Section 6 concludes. # Theory ## 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.[^1] 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}$. [^1]: Of course, while $y_t$ denotes daily new cases here, it could equally denote weekly sales of a new product, etc. None of the analysis here is dependent on the data frequency. ## 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 @Helske2017. 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 @Helske2017. 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.[^2] [^2]: These are not proper prediction or confidence intervals for $y_{T+\ell}$. The one-step-ahead predictive distribution of $\hat{y}_{T+\ell|T}$ (for $\ell=1$) is lognormal. This is not the case more than one step ahead due to the presence of the cumulative total in equation (2). 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 @Harvey2021a). 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[^3] [^3]: Note the change in sign of $\gamma$ as compared with HK. $$\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). @WallingaLipsitch2007 describe how the estimates of $g_{y,t}$ can be translated into estimates of the instantaneous reproduction number $R_{t}.$ @Harvey2021a 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} \noindent 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 @Harvey2021b 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. # Reinitialization The coronavirus pandemic was characterized by multiple waves punctuated by plateaus. At the beginning of a new wave the growth rate of daily cases, $g_{y,{t}}$, turns positive. The initial surge may be explosive to the point where the growth is *super exponential*. In this case, $\gamma _{t}$, the growth rate of $g_t$ (which is the growth rate of cumulative cases) can also turn positive, with no peak in prospect for $y_{t}.$ Such a phase can be expected to be transient, with $\gamma _{t}$ dropping back to zero (exponential growth in infection, accompanied by an upcoming peak in $y_{t}$), and then falling below zero (sub-exponential growth in infection). From the point-of-view of forecasting an epidemic, a peak must be in prospect even if it can only be expected some way into the future. There is thus a need for a solution to the problem of the estimated $\hat\gamma _{t\mid t}$ rising to positive values as it adapts to the upward surge in $y_{t}$, and remaining positive for any protracted period. This upward shift in $\hat\gamma _{t\mid t}$ can be averted by *reinitializing* the $\ln g_{t}$ series at the start of a new wave. This involves setting the cumulative total of cases $Y_t$ back to zero at, or around, the start of a new wave and setting $\gamma_t$ to zero so as to impose exponential growth. From the point-of-view of the relationship $g_{y_{t}} = g_{t} + \gamma _{t},$ the re-initialization effectively shifts *surplus* $\gamma_t$ emanating from super-exponential growth, into $\delta_t$ and therefore into $g_t$ (since $g_t = \exp \delta_t$). Note that, on the date of the re-initialization, $g_{y,t} = g_t$, since $\gamma_t=0$, and both $g_{y,t}$ and $g_t$ will be high because a new wave is taking off. ## 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. # Functionality of `tsgc` The main Gompertz class in `tsgc` is `SSModelDynamicGompertz`, which implements the models described in (2)-(4), with or without reinitialization. A model object is constructed with `SSModelDynamicGompertz$new(...)` or `tsgc::SSModelDynamicGompertz(...)` and its unknown parameters are estimated with the `estimate()` method, invoked as `estimate(model)`, which returns an object of the `FilterResults` class. This class wraps the underlying `KFAS` `KFS` object, together with additional methods for prediction, and an optional `idx_calendar` (see below) attached. The `SSModelDynamicGompertz` needs only a cumulative series `Y` as an input, supplied as an `idx_series` object -- an integer-position representation of a time series, consisting of a vector or matrix of values together with a `start` position. In our application, `Y` is the cumulative number of new coronavirus cases. The estimation window is specified directly on the model object via the `start` and `end` arguments, which are integer positions in `Y` rather than calendar dates. There is an option to specify the signal-to-noise ratio $q_{\zeta}$, rather than estimate it, and an option to specify the model to have a seasonal component through the `sea.period` argument. The period of a seasonal component is specified through the `sea.period` option. `idx_series` positions carry no calendar meaning of their own: they are plain integers, and the model, its estimation, and its forecasts are defined purely in terms of them. Real calendar time is reintroduced only as an optional, purely cosmetic layer for plotting, via an `idx_calendar` object passed through the `calendar` argument. An `idx_calendar` records an anchor position and the corresponding calendar date (or other real-world reference point), together with the step size and unit between consecutive positions; it is never consulted during estimation or filtering. Supplying a `calendar` therefore changes nothing about the fitted model -- only whether its plots are labelled with dates or with raw integer positions. Two helper functions, `xts_to_idx()` and `idx_to_pos()`, are provided to ease the transition from calendar-time data: the former converts an `xts`/`Date`-indexed series into an `idx_series` plus a matching `idx_calendar` in one step, and the latter translates a calendar date back into the integer position `idx_series` expects for `start`/`end`/`reinit.idx`. Both are demonstrated in the Illustration section below. Every calendar used in this vignette has a step size that is a single, uniform `amount`/`unit` pair; `idx_calendar` also supports compound and heterogeneous step sizes for series that need them, described in the Appendices "Compound calendar steps" and "Heterogeneous repeating cycles" below. The `SSModelDynamicGompertz` class allows the model to be estimated for a new wave without losing information from prior waves. It will accept the reinitialization position specified by `reinit.idx` (an integer position) or a `FilterResults` object from which it can extract the initial values. If the user wishes to reinitialize the model without using prior information (i.e. treat the new wave as an entirely separate epidemic), a reinitialization position can be specified through the `reinit.idx` option and `use.presample.info` can be set to `FALSE`. The `FilterResults` class contains prediction methods which can be applied to estimated dynamic Gompertz curve models (both reinitialized and non-reinitialized). `get_growth_y()` returns filtered or smoothed estimates of the growth rate of new cases ($g_t$), while `get_gy_ci()` returns the same with uncertainty bands. Forecasts of the incidence variable (new cases, $y_t$) can be obtained with the `predict_level()` call, and forecasts of all the states can be obtained with the `predict_all()` call. Several functions are available to generate plots of smoothed and filtered estimates and forecasts. `plot_log_forecast()` plots the forecast and realised values of $\ln(g_t)$. `plot_gy_components()` and `plot_gy_ci()` can be used to plot the smoothed or filtered growth rate, its components, and uncertainty bands, respectively. Forecasts of the incidence variable ($y_t$) and forecast intervals can be plotted using `plot_forecast()`, while `plot_holdout()` adds plots of prediction intervals and of realized outcomes over a holdout period to help evaluate forecast accuracy. `plot_compare_forecast()` overlays the forecasts from a list of fitted models -- for example a reinitialized and non-reinitialized fit, or a dynamic Gompertz model against a leading-indicator model -- on the same axes for visual comparison. Finally, the `reinitialise_dataframe()` function can be used to reinitialise a data series at a given `reinit.idx`. Every one of these plotting functions accepts an `axis` argument, of class `idx_axis_opts`, which controls how integer positions are translated to the x-axis and is described in full in the Appendix "Controlling the plot x-axis" below. Left at its default (`NULL`), the x-axis uses real calendar dates whenever a suitable `calendar` is attached to the model or results object, and plain integer positions otherwise. More details on how to use the methods and functions described are presented in the following section. ## 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. # Forecast accuracy and rolling-origin validation Forecast accuracy can be evaluated with a variety of scale-dependent and scale-free measures. In the notation of the package, the mean absolute percentage error is $$\mathrm{MAPE}=\frac{100}{n}\sum_{t=1}^n \left|\frac{y_t-\hat y_t}{y_t}\right|,$$ while the symmetric MAPE, as implemented in the package (scaled 0-100 rather than the more common 0-200 convention that includes a factor of 2 in the numerator), is defined as $$\mathrm{sMAPE}=\frac{100}{n}\sum_{t=1}^n \frac{|y_t-\hat y_t|}{|y_t|+|\hat y_t|},$$ and the mean absolute error and root mean squared error are given by $$\mathrm{MAE}=\frac{1}{n}\sum_{t=1}^n |y_t-\hat y_t|, \qquad \mathrm{RMSE}=\sqrt{\frac{1}{n}\sum_{t=1}^n (y_t-\hat y_t)^2}.$$ sMAPE is reported as the primary percentage-based comparison measure in the empirical illustrations above, since MAPE can become unstable when observed values are close to zero. The rolling-origin procedure used by `cross_val()` re-estimates the model at successive evaluation positions, holds back a fixed forecast horizon, and computes the chosen loss metric over that out-of-sample window, so that forecast accuracy is assessed under a design that matches the practical forecasting problem. # Illustration of the `tsgc` package In this section we provide a full working example of the `tsgc` package in `R` which implements the modeling framework for time series growth curves-based epidemic forecasting. `tsgc` comes with two example data sets relating to COVID-19: one for Gauteng province in South Africa (sourced from South Africa's official coronavirus online news and information [portal](https://sacoronavirus.co.za/)) and another for England (sourced from the official UK government [dashboard](https://ukhsa-dashboard.data.gov.uk/) for data and insights on coronavirus). In the example that follows, we use the data on confirmed cases in Gauteng. The data series is in cumulative form and ships with the package as an `xts` object with a `Date` index, as follows. ```{r load-gauteng-xts} data(gauteng, package = "tsgc") ``` The modelling code works with `idx_series` (integer-position data) rather than directly with `xts` objects, so that estimation and forecasting are defined purely in terms of position, independent of calendar frequency or gaps. Calendar time is reintroduced only for plotting, via an optional `idx_calendar` that records how integer positions map back to real dates. `xts_to_idx()` performs this conversion in one step for the common case of a series with no gaps, returning both the `idx_series` and a matching `idx_calendar`; `idx_to_pos()` is its companion for translating a calendar date into the integer position that `start`/`end`/`reinit.idx` expect. A minimal example of converting a `CSV`-derived data frame to this representation is given in the Appendix; the pattern below extends directly to that case. ```{r convert-gauteng} conv <- xts_to_idx(gauteng) gauteng_idx <- conv$series gauteng_cal <- conv$calendar gauteng_cal$anchor_name <- "first recorded case" ``` New COVID-19 cases reported for Gauteng province and their centered 7-day moving average presented in Figure 1 show a sequence of four waves over the period between 10 March 2020 and 5 January 2022. ```{r make-fig-1, fig.cap="Figure 1: New Cases and their centered 7-day moving average for Gauteng province in South Africa between 10 March 2020 and 5 January 2022."} mod_all <- SSModelDynamicGompertz(Y = gauteng_idx, calendar = gauteng_cal) plot(mod_all, series.name = "daily cases", MA_period = 7) ``` ## 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. ```{r set_options} 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. ```{r model_est_free_q} 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. ```{r model_est_fix_q} 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. ```{r model-diagnostics} summary(res) tsgc::print_model_diagnostics(res) ``` ## 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`. ```{r fig-2, fig.cap="Figure 2: Fourteen-day forecast of $ln(g_t)$ from 4 May 2021, for Gauteng province in South Africa."} tsgc::plot_log_forecast( res = res, Y = Y, n.ahead = n.forecasts, plt.start = tail(res$index, 1) - plt.length ) ``` From these results we can recover the forecasts of new cases from 4 May 2021, with their forecast bands. ```{r fig-3, fig.cap = "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."} tsgc::plot_forecast( res = res, n.ahead = n.forecasts, confidence.level = confidence.level, plt.start = tail(res$index, 1) - plt.length, series.name = "cases" ) ``` 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. ```{r fig-4, fig.cap = "Figure 4: Accuracy of the fourteen-day forecast of new cases from 20 April 2021 for Gauteng province in South Africa."} # 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 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. ```{r calc_rt} gen_int <- 4 r.t <- estimate_r0(res, gen_int, n.ahead = 7) r.t ``` ```{r fig-5, fig.cap = "Figure 5: Reproduction numbers for the 7-day period to 3 May 2021, for Gauteng province in South Africa."} plot_r0(res, gen_int = gen_int, n.ahead = 7, title = "Reproduction numbers for the 7-day period to 3 May 2021") ``` 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. ```{r write_results, eval=F} 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. ```{r change-est-date-end} 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. ```{r fig-6, fig.cap = "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."} # 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")) ``` The position for reinitialization is set accordingly, translating the reset date to an integer position via `idx_to_pos()`. ```{r set_reinit_date} 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. ```{r reinit_estim} 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. ```{r fig-7, fig.cap = "Figure 7: Forecast of $ln(g_t)$ after reinitialization."} 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.' ) ``` 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. ```{r fig-8, fig.cap = "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, Y = Y, n.ahead = 14, confidence.level = 0.68, series.name = "cases" ) ``` ```{r fig-9, fig.cap = "Figure 9: Forecast accuracy of the model with 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" ) ``` We can also compare the standard (no reinitialization) and reinitalized forecasts in a single plot: ```{r fig-10, fig.cap = "Figure 10: Comparison of standard (green) and reinitialized (blue) forecasts"} tsgc::plot_compare_forecast( results = list(res.reinit, res), n.ahead = 14, actual = Y ) ``` 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. ```{r li-setup} 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. ```{r li-plot, fig.cap = "Figure 11: Daily COVID cases and hospitalisations, England."} 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 ) ``` We then estimate the baseline leading-indicator model, with cases as the lead series (`LeadIndCol = 1`) and weekly seasonality. ```{r li-estimate} 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. ```{r li-diagnostics} summary(res.eng) tsgc::print_model_diagnostics(res.eng) ``` The fitted object `res.eng` can be used to generate the same forecast and holdout panels as for the dynamic Gompertz model. ```{r fig-li-log, fig.cap = "Figure 12: Forecast of log growth rate of hospitalisations, England."} 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)" ) ``` ```{r fig-li-forecast, fig.cap = "Figure 13: Forecast 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)" ) ``` ```{r fig-li-holdout, fig.cap = "Figure 14: Accuracy of the seven-day 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)" ) ``` ### 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. [^4]: Note that these weather variables are realised weather covering the forecast horizon being predicted, rather than an archived forecast available at the time of the forecast origin which could potentially induce look-ahead bias. This exercise is purely to illustrate the capabilities of the `tsgc` package. ```{r li-xpred} 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 ``` ```{r fig-li-x-forecast, fig.cap = "Figure 15: Forecast of hospitalisations with weather regressors (weather, oracle/realised), England."} 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" ) ``` 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. ```{r cv-setup} 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) ``` ```{r fig-uk-italy, fig.cap = "Figure 16: Daily COVID cases in UK and Italy."} 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 ) ``` 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. ```{r uk-gompertz} 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) ``` ```{r uk-lead} 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) ``` ```{r fig-uk-compare, fig.cap = "Figure 17: Comparison of dynamic Gompertz and leading-indicator forecasts, UK daily cases."} tsgc::plot_compare_forecast( list(res.uk.gomp, res.uk.lead), actual = Yuk, n.ahead = n.forecasts.uk ) ``` ```{r fig-uk-holdout, fig.cap = "Figure 18: Accuracy of the leading-indicator forecast, 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" ) ``` In this first-wave window, the leading-indicator model gives substantially smaller holdout errors than the UK-only dynamic Gompertz model. ### Selecting the lag by rolling cross-validation The baseline lag $k=14$ is motivated by the observed first-wave timing. To avoid selecting the lag purely by inspection, the package provides a rolling cross-validation utility, `cross_val()`, which evaluates candidate specifications over a rolling origin and reports a chosen accuracy criterion. Because `SSModelLeadingIndicator` uses a single alignment lag per fitted model, this should be read as model selection over alternative alignments, not as estimation of a distributed-lag response. The following cross-validation loop tests lags up to $k=21$. We must take care to set the estimation end date in the selection part to give us a sufficiently long window to estimate the models with longer lags in the cross-validation selection part. This estimation window ends later than the estimation window used for this data earlier in this vignette. ```{r cv-models} 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 ) } ``` `cross_val()` is run twice with disjoint sets of forecast origins, spaced by the forecast horizon itself so evaluation windows do not overlap: an earlier selection block used only to choose the lag, and a later reporting block used only to score the selected lag against the benchmarks. Using the same folds for both selection and reporting would overstate the selected model's accuracy. ```{r cv-run} 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 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 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 ``` For each candidate model, `cross_val()` evaluates forecasting performance over a rolling-origin scheme, re-estimating the model at each origin and recording the chosen accuracy criterion over the forecast horizon. `Naive_last_value` and `RW_growth` are naive benchmarks; the `Vanilla` rows are UK-only dynamic Gompertz benchmarks (fixed $q_\zeta = 0.005$ and the AR(1) extension respectively); the remaining rows are Italy-to-UK leading-indicator models at different lags. `cv_reporting` is the accuracy figure that should be quoted for the selected lag, since it comes from origins that played no role in choosing it; `cv_selection` should not itself be read as an accuracy claim about the winning lag. Note that this comparison still uses only a handful of origins from one epidemic wave and realised (not archived) values throughout, so it remains illustrative of the `cross_val()` API rather than a general validation of any one lag choice. ## 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. ```{r freq-quarterly} 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) ``` ```{r fig-freq-q, fig.cap = "Figure 19: Log forecasts of quarterly Switch sales, using Wii sales as a leading indicator."} tsgc::plot_log_forecast( res.q, Y = nintendo_idx, n.ahead = 8, title = "Log forecasts of Switch sales" ) ``` ```{r freq-monthly} 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) ``` ```{r fig-freq-m, fig.cap = "Figure 20: Forecasts of monthly AvaTrade downloads, using DEGIRO downloads as a leading indicator."} tsgc::plot_forecast( res.m, n.ahead = 4, title = "AvaTrade monthly downloads", series.name = "Downloads" ) ``` 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. # Conclusions The `tsgc` package is based on a dynamic Gompertz curve model for the log of the growth rate of cumulative cases in an epidemic, with seasonal terms that capture day-of-the-week effects. The estimation is carried out using `KFAS`, a package for state-space modeling in `R`. The Kalman filter is used to estimate the state vector at each time point. The filter is initialized using a diffuse prior for the initial state vector. We allow the signal-to-noise ratio to either be estimated or fixed at some value based on experience and judgment. Future observations are forecast using predictive recursions. Epidemics are often characterized by multiple waves. A natural problem in this context is that there is very little data pertinent to the new wave available in its initial stages. The package employs a reinitialization method using priors in a way that allows data from before the beginning of the new wave to be used in estimation. The same framework can be extended to a leading-indicator formulation for outcomes that are preceded by an observable correlate, and it can be applied to data measured at different frequencies once the relevant seasonal period and lag structure are chosen. The package is demonstrated using COVID-19 data from South Africa and England, but it can be used to model and forecast any time series variable where a growth curve-like trajectory is expected. Examples might include sales of a new product, innovation adoption or website traffic. The approach is phenomenological rather than mechanistic: it is designed for short-horizon forecasting and exploratory analysis rather than for causal identification or structural epidemiological inference. The forecast bands are state-implied summaries of the latent-state uncertainty and should not be interpreted as calibrated prediction or confidence intervals in the usual sense. In addition, the implied lead–lag relationship between a leading indicator and a target series can break down over time, so the model should be monitored and re-estimated as conditions change. # Acknowledgements {.unnumbered} We thank The Cambridge Centre for Health Leadership & Enterprise, University of Cambridge Judge Business School, and Public Health England/UK Health Security Agency for generous support. We are indebted to Thilo Klein and Stefan Scholtes for constructive comments. Andrew Harvey's work was carried out as part of the University of Cambridge Keynes Fund project \`Forecasting and Policy in Environmental Econometrics'. # References {.unnumbered} ::: {#refs} ::: # Appendix: Incorporating seasonal terms into the state space model When we add a seasonal term to the model, the observation equation in the dynamic Gompertz curve (2) becomes $$\ln g_{t}=\delta_{t}+\nu_t+\varepsilon_{t},\;\;\;\;\;\varepsilon_{t}\sim NID(0,\sigma_{\varepsilon }^{2}),\;\;\;\;\;t=s,...,T,$$ where $\nu_t$ is the seasonal component, $\delta_t$ remains defined by (3), and $\varepsilon_t$ remains the iid Normal disturbance. There are two options for specifying the evolution of the seasonal component. We can either use a trigonometric seasonal or a dummy variable seasonal. In our application, we use a trigonometric seasonal, although the two specifications are closely related. Indeed, @Proietti2000 shows that, under certain conditions, the two approaches are equivalent. In the trigonometric seasonal approach, the seasonal pattern is captured by a set of trigonometric terms at the seasonal frequencies $\lambda_j = \frac{2 \pi j}{s}$ for $j=1, \ldots s^*$, where $s^*=\frac{s}{2}$ if $s$, the periodicity of the seasonal effect, is even, and $s^*=\frac{s-1}{2}$ if $s$ is odd @Durbin2012. Our applications use daily data and we set $s=7$ to capture day-of-the-week effects. Letting $\nu_{j,t}$ be the effect of season $j$ at time $t$, the seasonal terms evolve according to \begin{align} \nu_{t} &= \sum_{j=1}^{s^{*}} \nu_{j,t}, & (13) \end{align} where \begin{align} \nu_{j,t} &= \nu_{j,t-1} \cos \lambda_j + \nu^*_{j,t-1} \sin \lambda_j + \omega_{j,t} & (14)\\ \nu^*_{j,t} &= -\nu_{j,t-1} \sin \lambda_j + \nu^*_{j,t-1} \cos \lambda_j + \omega^*_{j,t}, \;\;\;\;\; j = 1, \ldots, s^*, & (15) \end{align} and $\omega_{j,t}$ and $\omega^*_{j,t}$ are mutually independent, iid $N(0,\sigma^2_{\omega})$ variables. When reinitializing the model with seasonal terms, $P^r_1$ is a block-diagonal matrix based on $P_{r+1}$ which sets the covariances between $(\delta_t,\gamma_t)'$ and $(\nu_{1,t},\nu_{2,t},\ldots,\nu_{s^*,t})'$ to be zero. The covariance between $\delta_t$ and $\gamma_t$, as well as the covariances between $\nu_{1,t},\nu_{2,t},\ldots,\nu_{s^*,t}$, are permitted to be non-zero and come directly from $P_{r+1}$. # Appendix: How to import data from CSV file For models that include external regressors, forecasts are conditional on future values of those regressors. If the future regressor values originate in a CSV file with a `Date` column, the same `xts_to_idx()` conversion used elsewhere in this vignette applies: read the CSV, convert the date column, coerce to `xts`, and then convert to `idx_series` using the model's existing calendar so that the resulting positions line up correctly. The following example illustrates this using an inline text block for reproducibility; in practice the data would normally be read directly from a CSV file. The examples below assume a Gompertz model has already been fit with a weather regressor supplied via `xpred`, producing a `FilterResults` object called `res_weather`. This is fit explicitly here (using its own pinned estimation window, rather than the ambient `estimation.pos.end` reassigned above), so the CSV examples that follow are self-contained. ```{r fit-res-weather} data(gauteng_weather_2021, package = "tsgc") gauteng_weather_idx <- xts_to_idx( gauteng_weather_2021[, c(1, 3)], start.pos = idx_to_pos(cal, zoo::index(gauteng_weather_2021)[1]) )$series weather.pos.start <- estimation.pos.start weather.pos.end <- idx_to_pos(cal, "2021-05-03") gauteng_weather_est <- get_timeframe( gauteng_weather_idx, weather.pos.start, weather.pos.end ) model_weather <- SSModelDynamicGompertz( Y = Y, xpred = gauteng_weather_est, q = q, start = weather.pos.start, end = weather.pos.end, calendar = cal ) res_weather <- estimate(model_weather) ``` ```{r xpred-csv-example, eval = F} # This inline CSV block demonstrates the expected format of a file # containing future values of the external regressors. txt <- " Date,windspd_mtrs_p_sec,temperature_C 2021-05-04,2.03,14.61 2021-05-05,1.54,15.51 2021-05-06,1.94,16.42 2021-05-07,2.38,15.54 2021-05-08,2.57,14.18 2021-05-09,2.65,13.55 2021-05-10,2.19,14.84 2021-05-11,2.08,15.55 2021-05-12,2.07,15.97 2021-05-13,1.92,15.86 2021-05-14,1.85,15.79 2021-05-15,2.86,15.57 2021-05-16,3.46,16.42 2021-05-17,2.39,13.53 " # Read the demonstration CSV text into a data frame. gauteng_weather_future_csv <- read.csv( text = txt, stringsAsFactors = FALSE ) gauteng_weather_future_csv$Date <- as.Date(gauteng_weather_future_csv$Date) # Convert the data frame to xts first... gauteng_weather_future_xts <- xts::xts( gauteng_weather_future_csv[, -1], order.by = gauteng_weather_future_csv$Date ) # ...then to an idx_series, anchored on the fitted model's own calendar # so that the resulting positions line up with res_weather$calendar. gauteng_weather_future_idx <- xts_to_idx( gauteng_weather_future_xts, start.pos = idx_to_pos(res_weather$calendar, zoo::index(gauteng_weather_future_xts)[1]) )$series # Supply the future regressors directly on the fitted results object. res_weather$xpred.new <- gauteng_weather_future_idx ``` In practice, users would normally read the future regressor values directly from a CSV file rather than embedding the data in the code. The CSV file should have the same structure as shown above, with a `Date` column followed by one column for each external regressor. ```{r xpred-csv-file, eval = F} # Read future regressors from a CSV file. gauteng_weather_future_csv <- read.csv( "gauteng_weather_future.csv", stringsAsFactors = FALSE ) gauteng_weather_future_csv$Date <- as.Date(gauteng_weather_future_csv$Date) # Convert to xts, then to idx_series on the model's calendar. gauteng_weather_future_xts <- xts::xts( gauteng_weather_future_csv[, -1], order.by = gauteng_weather_future_csv$Date ) gauteng_weather_future_idx <- xts_to_idx( gauteng_weather_future_xts, start.pos = idx_to_pos(res_weather$calendar, zoo::index(gauteng_weather_future_xts)[1]) )$series # Supply the future regressors directly on the fitted results object. res_weather$xpred.new <- gauteng_weather_future_idx ``` Once the future regressors have been supplied, subsequent forecasts from the fitted model are conditional on these external regressor values. While this example is specifically for exogenous predictors supplied to a `SSModelDynamicGompertz` fit, the same `xts_to_idx()`/field-assignment pattern extends to `xpred_lead.new` and `xpred_targ.new` on a `SSModelLeadingIndicator` fit, as demonstrated for the England cases-hospitalisations example above. # Appendix: Controlling the plot x-axis with `idx_axis_opts()` Every plotting function in this section so far has left its `axis` argument unset, which defaults to `mode = "auto"`: real calendar dates if the calendar has `posixct = TRUE` (as in almost every example above), otherwise plain integer positions. `idx_axis_opts()` lets you override this per plot, and applies uniformly across every `plot.*`/`plot_*` function in the package. To keep the comparison clear, this section reuses a single plot and simply varies `axis`. ```{r axis-showcase-setup} # res/model/estimation.pos.end are re-used, generically-named objects # earlier in this vignette. res.axis.demo <- estimate( SSModelDynamicGompertz$new(Y = Y, q = q, sea.period = 7, start = estimation.pos.start, end = idx_to_pos(cal, "2021-05-03"), calendar = cal) ) showcase_axis_plot <- function(axis = NULL) { tsgc::plot_forecast( res = res.axis.demo, n.ahead = n.forecasts, confidence.level = confidence.level, plt.start = tail(res.axis.demo$index, 1) - plt.length, series.name = "cases", axis = axis ) } ``` ## `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. ```{r fig-axis-date, fig.cap = "Real calendar dates (mode = \"date\")."} showcase_axis_plot(idx_axis_opts(mode = "date")) ``` ## `mode = "position"` The raw integer `idx_series` position, ignoring the calendar entirely. ```{r fig-axis-position, fig.cap = "Raw integer positions (mode = \"position\")."} showcase_axis_plot(idx_axis_opts(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. ```{r fig-axis-steps, fig.cap = "Steps from the anchor (mode = \"steps\")."} showcase_axis_plot(idx_axis_opts(mode = "steps")) ``` ## `mode = "time_since"` The pattern-weighted calendar offset from the anchor, expressed in `calendar$unit`s 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. ```{r fig-axis-time-since, fig.cap = "Time since the anchor, in calendar units (mode = \"time_since\")."} showcase_axis_plot(idx_axis_opts(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. ```{r fig-axis-infobox, fig.cap = "Steps from anchor, with an info box summarising the calendar."} showcase_axis_plot(idx_axis_opts(mode = "steps", info_box = TRUE)) ``` If the info box's pattern were long (e.g. a business-day calendar with weekends skipped), `pattern_n` keeps the caption readable: ```{r fig-axis-infobox-truncated, fig.cap = "Info box with pattern truncated to its first 3 values via pattern_n."} showcase_axis_plot(idx_axis_opts(mode = "steps", info_box = TRUE, pattern_n = 3)) ``` # Appendix: Compound calendar steps with `idx_step()` and `idx_calendar_step()` Every calendar used elsewhere in this vignette is built with the primary `idx_calendar()` constructor, whose step size is a single `amount`/`unit` pair (e.g. `amount = 1, unit = "days"`, or `amount = 1, unit = "quarters"`). This covers the common case where the step between consecutive `idx_series` positions is a whole multiple of one calendar unit. Some series, however, are genuinely spaced by a combination of units that a single `amount`/`unit` pair cannot express - for example, a reporting system that timestamps quarterly figures a few seconds after the quarter boundary, so that the step is "one quarter plus three seconds". `idx_step()` builds this kind of compound step, and `idx_calendar_step()` is the `idx_calendar` constructor that takes one. ```{r appendix-idx-step} qtr_step <- idx_step(quarters = 1, seconds = 3) qtr_step cal_step <- idx_calendar_step( anchor = as.POSIXct("2024-01-01 00:00:03", tz = "UTC"), step = qtr_step, posixct = TRUE ) idx_to_date(cal_step, 1:4) ``` `idx_to_pos()` inverts this exactly as it does for the primary constructor. A compound step has no single `amount` to divide by, so the position is located by a binary search on the number of steps from the anchor, as it is for `idx_calendar()` calendars with numeric patterns: ```{r appendix-idx-step-inverse} idx_to_pos(cal_step, idx_to_date(cal_step, 4)) ``` Any component of an `idx_step` (`years`, `quarters`, `months`, `weeks`, `days`, `hours`, `minutes`, `seconds`) left at its default of `0` is simply omitted from the step, so `idx_step(days = 1)` is equivalent to the `amount = 1, unit = "days"` pairs used throughout the main text. `idx_step_add()` is the lower-level function that applies a given number of whole steps to an anchor date directly, and is what an `idx_calendar_step()`-built calendar uses internally to convert positions to dates; it is not usually called directly by users of the package; it is available for cases where a step needs to be applied outside of an `idx_calendar` altogether. # Appendix: Heterogeneous repeating cycles with `multi_step_pattern()` and `idx_calendar_multi_step()` The `pattern` argument of `idx_calendar()`/`idx_calendar_step()` handles a repeating cycle of *different multiples* of one step - e.g. a business-day calendar, where four single-day steps (Mon-Thu) are followed by a three-day step over the weekend. That still assumes every step in the cycle is the same *kind* of step, just scaled differently. Some data instead cycle through genuinely different kinds of steps - for example, "three days, three days, then one month, repeating". `multi_step_pattern()` describes this kind of cycle as a sequence of `idx_step` objects, and `idx_calendar_multi_step()` is the `idx_calendar` constructor that uses it. ```{r appendix-multi-step} cal_multi <- idx_calendar_multi_step( anchor = as.Date("2024-01-01"), multi_step = multi_step_pattern( idx_step(days = 3), idx_step(days = 3), idx_step(months = 1) ), posixct = TRUE ) idx_to_date(cal_multi, 1:6) ``` Because the slots in a `multi_step_pattern` can mix calendar-relative (e.g. months) and fixed-duration (e.g. days) components, there is no closed-form way to collapse the whole cycle into a single numeric offset. `idx_to_date()` instead walks one slot at a time from the anchor position, and `idx_to_pos()`'s inverse walks the same cycle in reverse: ```{r appendix-multi-step-inverse} idx_to_pos(cal_multi, "2024-02-07") ``` This makes `idx_calendar_multi_step()` calendars more expensive to convert than the other two constructors - each position takes time proportional to its distance from the anchor, rather than constant time - so it is best reserved for genuinely irregular step sequences that the `pattern` argument of the other two constructors cannot express. The `"time_since"` axis mode described in the previous appendix is also unavailable for this calendar: a `multi_step_pattern` has no single `amount` to express an offset in, so `mode = "steps"` or `mode = "date"` should be used instead. # Appendix: Non-calendar `idx_calendar` anchors and `idx_offset_to_pos()` Not every `idx_series` needs a genuine calendar interpretation. An `idx_calendar`'s `anchor` can be any single reference point, not just a `Date`/`POSIXct` - for example, a plain numeric count of picoseconds since the start of an experiment. `idx_offset_to_pos()` is the inverse operation for this case: given a value already expressed in the anchor's own units, it returns the integer `idx_series` position that corresponds to it - the counterpart of `idx_to_pos()`, which instead expects a calendar date. ```{r appendix-offset-to-pos} cal_ps <- idx_calendar(anchor = 0, anchor_pos = 1L, amount = 2.5, unit = "picoseconds") idx_offset_to_pos(cal_ps, 12.5) ``` None of the examples elsewhere in this vignette need `idx_offset_to_pos()`, since every series used above has a genuine `Date`/`POSIXct` calendar and so uses `idx_to_pos()` instead; it is included here for completeness, for applications - e.g. raw instrumentation data with no calendar meaning - where `idx_series` positions are still worth translating to and from the anchor's own numeric scale.