Movatterモバイル変換

FilterResults

Description

Class for estimated Dynamic Gompertz Curve model and containsmethods to extract smoothed/filtered estimates of the states, the level ofthe incidence variabley, and forecasts ofy.

Methods

get_growth_y(smoothed = FALSE, return.components = FALSE)

Returns the growth rate of the incidence (y) of the cumulatedvariable (Y). Computed as

g_t = \exp\{\delta_t\}+\gamma_t.

Parameters

• smoothed Logical value indicating whether to use thesmoothed estimates of\delta and\gamma to compute thegrowth rate (TRUE), or the contemporaneous filtered estimates(FALSE). Default isFALSE.

• return.components Logical value indicating whether toreturn the estimates of\delta and\gamma as well asthe estimates of the growth rate, or just the growth rate. Default isFALSE.

Return Value

xts object containingsmoothed/filtered growth rates and components (\delta and\gamma), where applicable.

get_gy_ci(smoothed = FALSE, confidence_level = 0.68)

Returns the growth rate of the incidence (y) of the cumulatedvariable (Y). Computed as

g_t = \exp\{\delta_t\}+\gamma_t.

Parameters

• smoothed Logical value indicating whether to use thesmoothed estimates of\delta and\gamma to compute thegrowth rate (TRUE), or the contemporaneous filtered estimates(FALSE). Default isFALSE.

• confidence_level Confidence level for the confidenceinterval. Default is0.68, which is one standard deviation fora normally distributed random variable.

Return Value

xts object containing smoothed/filteredgrowth rates and upper and lower bounds for the confidence intervals.

predict_all(n.ahead, sea.on = FALSE, return.all = FALSE)

Returns forecasts of the incidence variabley, the state variablesand the conditional covariance matrixfor the states.

Parameters

• n.ahead The number of forecasts you wish to create fromthe end of your sample period.

• sea.on Logical value indicating whether seasonalcomponents should be included in thestate-space model or not. Default isTRUE.

• return.all Logical value indicating whether to returnall filtered estimates and forecasts(TRUE) or only the forecasts (FALSE). Default isFALSE.

Return Value

xts object containing the forecast(and filtered, where applicable) levelofy (y.hat),\delta (level.t.t),\gamma (slope.t.t), vector of states including theseasonals where applicable (a.t.t) and covariance matrix of allstates including seasonals where applicable (P.t.t).

predict_level( y.cum, n.ahead, confidence_level, sea.on = FALSE, return.diff = FALSE)

Forecast the cumulated variable or the incidence of it. This function returnsthe forecast of the cumulated variableY, or the forecast of the incidence of the cumulated variable,y. Forexample, in the case of an epidemic,y might be daily new cases ofthe disease andY the cumulative number of recorded infections.

Parameters

• y.cum The cumulated variable.

• n.ahead The number of periods ahead you wish to forecast fromthe end of the estimation window.

• confidence_level The confidence level for the log growthrate that should be used to computethe forecast intervals ofy.

• return.diff Logical value indicating whether to return the cumulated variable,Y, or the incidence of it,y (i.e., the first difference of the cumulated variable). Default isFALSE.

Return Value

xts object containing the pointforecasts and upper and lower bounds ofthe forecast interval.

print_estimation_results()

Prints a table of estimated parameters in a format ready to paste intoLaTeX.

References

Harvey, A. C. and Kattuman, P. (2021). A Farewell to R:Time Series Models for Tracking andForecasting Epidemics, Journal of the Royal Society Interface, vol 18(182):20210179

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Print estimation resultsres$print_estimation_results()# Forecast 7 days ahead from the end of the estimation windowres$predict_level(y.cum = gauteng[idx.est], n.ahead = 7,  confidence_level = 0.68)# Forecast 7 days ahead from the model and return filtered statesres$predict_all(n.ahead = 7, return.all = TRUE)# Return the filtered growth rate and its componentsres$get_growth_y(return.components = TRUE)# Return smoothed growth rate of incidence variable and its confidence# intervalres$get_gy_ci(smoothed = TRUE, confidence_level = 0.68)

Base class for estimating time-series growth curve models. ClassesSSModelDynamicGompertz andSSModelDynGompertzReinit refer backto this base class.

Description

Base class for estimating time-series growth curve models. ClassesSSModelDynamicGompertz andSSModelDynGompertzReinit refer backto this base class.

Methods

estimate(sea.type = "trigonometric", sea.period = 7)

Estimates the dynamic Gompertz curve model when applied to an object ofclassSSModelDynamicGompertz orSSModelDynGompertzReinit.

Parameters

• sea.type Seasonal type. Options are'trigonometric' and'none'.'trigonometric' willyield a model with a trigonometric seasonal component and'none' will yield a model with no seasonal component.

• sea.period The period of seasonality. For aday-of-the-week effect with daily data, this would be 7. Not requiredifsea.type = 'none'.

Return Value

An object of classFilterResultscontaining the result output for the estimated dynamic Gompertz curvemodel.

get_dynamic_gompertz_model( y, q = NULL, sea.type = "trigonometric", sea.period = 7, a1 = NULL, P1 = NULL, Q = NULL, H = NULL)

Returns dynamic Gompertz curve model.

Parameters

• y The cumulated variable

• q The signal-to-noise ratio (ratio of slope to irregularvariance). Defaults to'NULL', in which case no signal-to-noiseratio will be imposed. Instead, it will be estimated.

• sea.type Seasonal type. Options are'trigonometric'and'none'.'trigonometric' will yield a model with atrigonometric seasonal component and'none' will yield a modelwith no seasonal component.

• sea.period The period of seasonality. For a day-of-the-weekeffect with daily data, this would be 7. Not required ifsea.type = 'none'.

• a1 Optional parameter specifying the prior mean of thestates. Defaults to'NULL'. Leave as'NULL' for a diffuseprior (no prior information). If a proper prior is to be specified, botha1 andP1 must be given.

• P1 Optional parameter specifying the prior mean of thestates. Defaults to'NULL'. Leave as'NULL' for a diffuseprior (no prior information). If a proper prior is to be specified,botha1 andP1 must be given.

• Q Optional parameter specifying the state error varianceswhere these are to be imposed rather than estimated. Defaults to'NULL' which will see the variances estimated.

• H Optional parameter specifying the irregular variancewhere this is to be imposed rather than estimated. Defaults to'NULL' which will see the variance estimated.

Description

The dynamic Gompertz with an integrated random walk (IRW) trend is

\ln g_{t}=\delta_{t}+\varepsilon_{t}, \;\;\;\; \varepsilon_{t}\sim NID(0,\sigma_{\varepsilon }^{2}), \;\;\;\; t=2,...,T,

whereY_t is the cumulated variable,y_t = \Delta Y_t,\ln g_{t}=\ln y_{t}-\ln Y_{t-1} and

\delta_{t} =\delta_{t-1}+\gamma_{t-1},

\gamma_{t} =\gamma_{t-1}+\zeta_{t}, \;\;\;\; \zeta_{t}\sim NID(0,\sigma_{\zeta }^{2}),

where the observation disturbances\varepsilon_{t} and slopedisturbances\zeta_{t}, are iid Normal and mutually independent.Note that, the larger the signal-to-noise ratio,q_{\zeta }=\sigma_{\zeta }^{2}/\sigma_{\varepsilon }^{2},the faster the slope changes in response to new observations. Conversely,a lower signal-to-noise ratio induces smoothness.

For the model without seasonal terms (sea.type = 'none') the arepriors are

\begin{pmatrix} \delta_1 \ \gamma_1 \end{pmatrix} \sim N(a_1,P_1)

.The diffuse prior hasP_1 = \kappa I_{2\times 2} with\kappa \to \infty. Implementation of the diffuse prior is handledby the packageKFAS (Helske, 2017). Where the model has a seasonalcomponent (sea.type = 'trigonometric'), the vector of prior meansa_1 and the prior covariance matrixP_1 are extendedaccordingly.

See the vignette for details of the variance matrixQ.H = \sigma^2_{\varepsilon}.

update(pars, model, q, sea.type)

Update method for Kalman filter to implement the dynamic Gompertz curvemodel.A maximum of 3 parameters are used to set the observation noise(1 parameter), the transition equation slope and seasonal noise. If q (signalto noise ratio) is not null then the slope noise is set using thisratio.

Parameters

• pars Vector of parameters.

• modelKFS model object.

• q The signal-to-noise ratio (ratio of slope to irregularvariance).

• sea.type Seasonal type. Options are'trigonometric' and'none'.

Return Value

KFS model object.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-06")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()

Class for re-initialised dynamic Gompertz curve model

Description

This class allows the implementation of the reinitialisationprocedure described in the vignette and summarised below.Lett=r denote the re-initialization date andr_0 denote thedate at which the cumulative series is set to 0. As the growth rate ofcumulative cases is defined asg_t\equiv \frac{y_t}{Y_{t-1}}, we have:

\ln g_t = \ln y_t - \ln Y_{t-1} \;\;\;\; t=1, \ldots, r

\ln g_t^r = \ln y_t - \ln Y_{t-1}^r \;\;\;\; t=r+1, \ldots, T

Y_{t}^{r}=Y_{t-1}^{r}+y_{t} \;\;\;\; t=r,\ldots,T

whereY_{t}^{r} is the cumulative cases after re-initialization. Wechoose to set the cumulative cases to zero atr_0=r-1, Y_{r-1}^{r}=0,such that the growth rate of cumulative cases is available fromt=r+1onwards.We reinitialise the model by specifying the prior distribution for theinitial states appropriately. See the vignette for details.

Methods

• new(Y, q = NULL, reinit.date=NULL, original.results=NULL,use.presample.info=TRUE) Create an instance of theSSModelDynGompertzReinit class.

  Parameters

  • Y The cumulated variable.

  • q The signal-to-noise ratio (ratio of slope to irregularvariance). Defaults to'NULL', in which case no signal-to-noiseratio will be imposed. Instead, it will be estimated.

  • reinit.date The reinitialisation dater. Should bespecified as an object of class"Date". Must be specified.

  • original.results Rather than re-estimating the model upto thereinit.date, aFilterResults class object can bespecified here and the parameters for the reinitialisation will be takenfrom this object. Default isNULL. This parameter is optional.

  • use.presample.info Logical value denoting whether ornot to use information from before the reinitialisation date in thereinitialisation procedure. Default isTRUE. IfFALSE, themodel is estimated from scratch from the reinitialisation date and noattempt to use information from before the reinitialisation date is made.

• get_model(y, q=NULL, sea.type = NULL, sea.period) Retrievesthe model object, which is a dynamic Gompertz curve model reinitialised atself$reinit.date.

  Parameters

  • y The cumulated variable.

  • q The signal-to-noise ratio (ratio of slope to irregularvariance). Defaults to'NULL', in which case no signal-to-noise ratiowill be imposed. Instead, it will be estimated.

  • sea.type Seasonal type. Options are'trigonometric' and'none'.'trigonometric' will yield a model with atrigonometric seasonal component and'none' will yield a model withno seasonal component.

  • sea.period The period of seasonality. For a day-of-the-weekeffect with daily data, this would be 7. Not required ifsea.type = 'none'.

  Return Value

  KFS model object.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2021-05-20")# Specify a modelmodel.reinit <- SSModelDynGompertzReinit$new(Y = gauteng[idx.est], q = 0.005,  reinit.date = as.Date("2021-04-29"))# Estimate a specified modelres.reinit <- model.reinit$estimate()## Alternatively, we could feed in a prior results object rather than a## reinitialisation date. The results are identical to the above.# Specify initial modelidx.orig <- zoo::index(gauteng) <= as.Date("2021-04-29")model.orig <- SSModelDynamicGompertz$new(Y = gauteng[idx.orig], q = 0.005)res.orig <- model.orig$estimate()# Estimate a specified modelmodel.reinit2 <- SSModelDynGompertzReinit$new(Y = gauteng[idx.est],q = 0.005, reinit.date = as.Date("2021-04-29"), original.results = res.orig)res.reinit2 <- model.reinit2$estimate()

Class for dynamic Gompertz curve state space model object.

Description

Class for dynamic Gompertz curve state space model object.

Methods

get_model(y, q = NULL, sea.type = 'trigonometric', sea.period = 7)Retrieves the model object.

Parameters

• y The cumulated variable.

• q The signal-to-noise ratio (ratio of slope to irregularvariance). Defaults to'NULL', in which case no signal-to-noise ratiowill be imposed. Instead, it will be estimated.

• sea.type Seasonal type. Options are'trigonometric' and'none'.'trigonometric' will yield a model with a trigonometricseasonal component and'none' will yield a model with no seasonalcomponent.

• sea.period The period of seasonality. For a day-of-the-weekeffect with daily data, this would be 7. Not required ifsea.type = 'none'.

Return Value

KFS model object.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-06")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()

Return index and value of maximum

Description

Something similar to Python's argmax.

Usage

argmax(x, decreasing = TRUE)

Arguments

Value

The maximum value and its index.

Examples

library(tsgc)data(gauteng,package="tsgc")argmax(gauteng)

Compute log growth rate of cumulated dataset

Description

Helper method to compute the log growth rates of cumulatedvariables. It will compute the log cumulative growth rate for each column inthe data frame.

Usage

df2ldl(dt)

Arguments

Value

A data frame of log growth rates of the cumulated variable which hasbeen inputted via the parameterdt.

Examples

library(tsgc)data(gauteng,package="tsgc")df2ldl(gauteng)

Cumulative cases of Covid-19 in England.

Description

Cumulative cases of Covid-19 in England.

Usage

data(england)

Format

An object of class"xts";

Cases

Cumulative cases of Covid-19

References

Downloaded from https://ukhsa-dashboard.data.gov.uk/topics/covid-19

Examples

data(england)# plot daily casesplot(diff(england))

Returns forecast of number of periods until peak given estimatedstate variables\delta and\gamma.

Description

Since Harvey and Kattuman (2021) show that

g_{y,t+\ell|T} = \exp\{\delta_{T|T}+\ell \gamma_{T|T}\}+\gamma_{T|T},

we can compute the\ell for whichg_{y,t}=0 and then will fallbelow zero. This\ell is given by

\ell = \frac{\ln(-\gamma_{T|T})-\delta_{T|T}}{\gamma_{T|T}}.

This ispredicated on\gamma_{T|T}<0, else there is super-exponential growth anno peak in sight. Of course, it only makes sense to investigate an upcomingpeak forg_{y,T|T}>0 (when cases are growing).

Usage

forecast.peak(delta, gamma)

Arguments

Value

Forecast of number of periods until peak.

Examples

# Forecasts the peak of an epidemic with gamma < 0 so that a peak is in# sight.forecast.peak(-2.87,-0.045)# Does not return a result (returns an error as gamma > 0)try(forecast.peak(-2.87,0.045), silent=TRUE)

Returns forecast of number of periods until peak givenKFAS::KFS output.

Description

Since Harvey and Kattuman (2021) show that

g_{y,t+\ell|T}= \exp\{\delta_{T|T}+\ell \gamma_{T|T}\}+\gamma_{T|T},

we can compute the\ell for whichg_{y,t}=0 and then will fall below zero. This\ell is given by

\ell = \frac{\ln(-\gamma_{T|T})-\delta_{T|T}}{\gamma_{T|T}}.

This ispredicated on\gamma_{T|T}<0, else there is super-exponential growthand no peak in sight. Of course, it only makes sense to investigate anupcoming peak forg_{y,T|T}>0 (when cases are growing). The estimatesof\delta_{T|T} and\gamma_{T|T} are extracted from theKFS object passed to the function.

Usage

forecast_peak(kfs_out)

Arguments

Value

Forecast of number of periods until peak.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-06")res <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)$estimate()forecast_peak(res$output)

Cumulative cases of Covid-19 in the South African province of Gauteng.

Description

Cumulative cases of Covid-19 in the South African province of Gauteng.

Usage

data(gauteng)

Format

An object of class"xts";

Cases

Cumulative cases of Covid-19 from 10th March 2020

References

Downloaded from https://sacoronavirus.co.za/

Examples

data(gauteng)# plot daily casesplot(diff(gauteng))

Plots forecast and realised values of the log cumulative growth rate

Description

Plots actual and filtered values of the log cumulative growthrate (\ln(g_t)) in the estimation sample and the forecast and realisedlog cumulative growth rate out of the estimation sample.

Usage

plot_forecast(  res,  y.eval,  n.ahead = 14,  plt.start.date = NULL,  title = "",  caption = "")

Arguments

Value

Aggplot2 plot.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")idx.eval <- (zoo::index(gauteng) >= as.Date("2020-07-20")) &     zoo::index(gauteng) <= as.Date("2020-07-27")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Plot forecast and realised log growth rate of cumulative casesplot_forecast(res, y.eval = df2ldl(gauteng[idx.eval]), n.ahead = 7,  title = "Forecast ln(g)", plt.start.date = as.Date("2020-07-13"))

Plots the growth rates and slope of the log cumulative growth rate

Description

Plots the smoothed/filtered growth rate of the difference in thecumulated variable (g_y) and the associated confidence intervals.

Usage

plot_gy_ci(  res,  plt.start.date = NULL,  smoothed = FALSE,  title = NULL,  series.name = NULL,  pad.right = NULL)

Arguments

Value

Aggplot2 plot.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Plot filtered gy, g and gammaplot_gy_ci(res, plt.start.date = as.Date("2020-07-13"))

Plots the growth rates and slope of the log cumulative growth rate

Description

Plots the smoothed/filtered growth rate of the difference inthe cumulated variable (g_y), the smoothed/filtered growth rate of thethe cumulated variable (g), and the smoothed/filtered slope of\ln(g),\gamma.Following Harvey and Kattuman (2021), we computeg_{y,t} as

g_{y,t} = \exp(\delta_t) + \gamma_t.

Usage

plot_gy_components(res, plt.start.date = NULL, smoothed = FALSE, title = NULL)

Arguments

Value

Aggplot2 plot.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Plot filtered gy, g and gammaplot_gy_components(res, plt.start.date = as.Date("2020-07-06"))

Plots the forecast of new cases (the difference of the cumulatedvariable) over a holdout sample.

Description

Plots actual values of the difference in the cumulated variable,the forecasts of the cumulated variable (both including and excluding theseasonal component, where a seasonal is specified) and forecast intervalsaround the forecasts, plus the actual outcomes from the holdout sample. Theforecast intervals are based on the prediction intervals for\ln(g_t).Also reports the mean absolute percentage prediction error over the holdoutsample.

Usage

plot_holdout(  res,  Y,  Y.eval,  confidence.level = 0.68,  date_format = "%Y-%m-%d",  series.name = NULL,  title = NULL,  caption = NULL)

Arguments

Value

Aggplot2 plot.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")idx.eval <- (zoo::index(gauteng) >= as.Date("2020-07-20")) &     zoo::index(gauteng) <= as.Date("2020-07-27")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Plot forecasts and outcomes over evaluation periodplot_holdout(res = res, Y = gauteng[idx.est], Y.eval = gauteng[idx.eval])

Plots the forecast of new cases (the difference of the cumulatedvariable)

Description

Plots actual values of the difference in the cumulated variable,the forecasts of the cumulated variable (both including and excluding theseasonal component, where a seasonal is specified) and forecast intervalsaround the forecasts. The forecast intervals are based on the predictionintervals for\ln(g_t).

Usage

plot_new_cases(  res,  Y,  n.ahead,  confidence.level = 0.68,  date_format = "%Y-%m-%d",  title = NULL,  plt.start.date = NULL)

Arguments

Value

Aggplot2 plot.

Examples

library(tsgc)data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-20")# Specify a modelmodel <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)# Estimate a specified modelres <- model$estimate()# Plot forecast of new cases 7 days aheadplot_new_cases(res, Y = gauteng[idx.est], n.ahead = 7,confidence.level = 0.68, date_format = "%Y-%m-%d",title = "Forecast new cases", plt.start.date = as.Date("2020-07-13"))

Reinitialise a data frame by subtracting thereinit.date row fromall columns

Description

Reinitialise a data frame by subtracting thereinit.date row fromall columns

Usage

reinitialise_dataframe(dt, reinit.date)

Arguments

Value

The reinitialised data frame

Examples

library(tsgc)data(gauteng,package="tsgc")reinitialise_dataframe(gauteng,as.Date("2021-01-01"))

Write a selection of relevant results to disc

Description

Function writes the following results to csv files which getsaved in the location specified inres.dir: forecast new cases orincidence variable,y; the filtered level and slope of\ln g,\delta and\gamma; filtered estimates ofg_y and theconfidence intervals for these estimates.

Usage

write_results(res, res.dir, Y, n.ahead, confidence.level)

Arguments

Value

A number of csv files saved in the directory specified inres.dir.

Examples

# Not run as do not wish to save to local disc when compiling documentation.# Below will run if copied and pasted into console.library(tsgc)library(here)res.dir <- tempdir()data(gauteng,package="tsgc")idx.est <- zoo::index(gauteng) <= as.Date("2020-07-06")res <- SSModelDynamicGompertz$new(Y = gauteng[idx.est], q = 0.005)$estimate()tsgc::write_results(res=res, res.dir = res.dir, Y = gauteng[idx.est], n.ahead = 14,confidence.level = 0.68)