`"hhh4"`

Count Time Series`hhh4_simulate.Rd`

Simulates a multivariate time series of counts based on the Poisson/Negative Binomial model as described in Paul and Held (2011).

```
# S3 method for hhh4
simulate(object, nsim = 1, seed = NULL, y.start = NULL,
subset = 1:nrow(object$stsObj), coefs = coef(object),
components = c("ar","ne","end"), simplify = nsim>1, ...)
```

- object
an object of class

`"hhh4"`

.- nsim
number of time series to simulate. Defaults to

`1`

.- seed
an object specifying how the random number generator should be initialized for simulation (via

`set.seed`

). The initial state will also be stored as an attribute`"seed"`

of the result. The original state of the`.Random.seed`

will be restored at the end of the simulation. By default (`NULL`

), neither initialization nor recovery will be done. This behaviour is copied from the`simulate.lm`

method.- y.start
vector or matrix (with

`ncol(object$stsObj)`

columns) with starting counts for the epidemic components. If`NULL`

, the observed means in the respective units of the data in`object`

during`subset`

are used.- subset
time period in which to simulate data. Defaults to (and cannot exceed) the whole period defined by the underlying

`"sts"`

object.- coefs
coefficients used for simulation from the model in

`object`

. Default is to use the fitted parameters. Note that the`coefs`

-vector must be in the same order and scaling as`coef(object)`

, which especially means`reparamPsi = TRUE`

(as per default when using the`coef`

-method to extract the parameters). The overdispersion parameter in`coefs`

is the inverse of the dispersion parameter`size`

in`rnbinom`

.- components
character vector indicating which components of the fitted model

`object`

should be active during simulation. For instance, a simulation with`components="end"`

is solely based on the fitted endemic mean.- simplify
logical indicating if only the simulated counts (

`TRUE`

) or the full`"sts"`

object (`FALSE`

) should be returned for every replicate. By default a full`"sts"`

object is returned iff`nsim=1`

.- ...
unused (argument of the generic).

Simulates data from a Poisson or a Negative Binomial model
with mean
$$\mu_{it} = \lambda_{it} y_{i,t-1} +
\phi_{it} \sum_{j \neq i} w_{ji} y_{j,t-1} +
\nu_{it}$$
where
\(\lambda_{it}>0\), \(\phi_{it}>0\), and \(\nu_{it}>0\) are
parameters which are modelled parametrically.
The function uses the model and parameter estimates of the fitted
`object`

to simulate the time series.

With the argument `coefs`

it is possible to simulate from
the model as specified in `object`

, but with different
parameter values.

If `simplify=FALSE`

: an object of class
`"sts"`

(`nsim = 1`

) or a list of those
(`nsim > 1`

).

If `simplify=TRUE`

: an object of class
`"hhh4sims"`

, which is an array of dimension
`c(length(subset), ncol(object$stsObj), nsim)`

, where the third
dimension is dropped if `nsim=1`

(yielding a matrix).
The originally observed counts during the simulation period,
`object$stsObj[subset,]`

, are attached for reference
(used by the `plot`

-methods) as an attribute `"stsObserved"`

,
and the initial condition `y.start`

as attribute `"initial"`

.

Paul, M. and Held, L. (2011) Predictive assessment of a non-linear
random effects model for multivariate time series of infectious
disease counts. Statistics in Medicine, **30**, 1118--1136

Michaela Paul and Sebastian Meyer

```
data(influMen)
# convert to sts class and extract meningococcal disease time series
meningo <- disProg2sts(influMen)[,2]
# fit model
fit <- hhh4(meningo, control = list(ar = list(f = ~ 1),
end = list(f = addSeason2formula(S = 1, period = 52)),
family = "NegBin1"))
plot(fit)
# simulate from model
simData <- simulate(fit, seed=1234)
# plot simulated data
plot(simData, main = "simulated data", xaxis.labelFormat=NULL)
# consider a Poisson instead of a NegBin model
coefs <- coef(fit)
coefs["overdisp"] <- 0
simData2 <- simulate(fit, seed=123, coefs = coefs)
plot(simData2, main = "simulated data: Poisson model", xaxis.labelFormat = NULL)
# consider a model with higher autoregressive parameter
coefs <- coef(fit)
coefs[1] <- log(0.5)
simData3 <- simulate(fit, seed=321, coefs = coefs)
plot(simData3, main = "simulated data: lambda = 0.5", xaxis.labelFormat = NULL)
```