Package {DynCount}

DynCount: Bayesian Dynamic Models for Poisson and Binomial Time Series

Description

The DynCount package fits GMRF state-space models for non-Gaussian time series. A latent trajectory z_t follows flexible dynamics – a first-order random walk or a stationary AR(1) process – and the observations are linked to it through either a Poisson (log link) or a binomial (logit link) observation model. An optional known offset may be added to the observation model's linear predictor (a log-exposure for the Poisson mean, a logit shift for the binomial probability); it is a fixed, user-supplied input and is not part of the latent process z_t.

Details

The package implements and extends the methodology of Zens and Bijak (2026, doi:10.1214/26-AOAS2171), adding stationary AR(1) dynamics, an optional drift/intercept, several innovation structures (Gaussian, Student-t, finite scale mixture, stochastic volatility), time-constant zero inflation, and fully propagated forecasting during MCMC.

Main entry points

fit_dynamic_model()

Fit a Poisson or binomial dynamic model (random walk or stationary AR(1)), optionally forecasting during sampling via horizon.

forecast() / predict()

Extract the forecast draws produced during sampling, or the in-sample fitted values.

summary()

Posterior summaries of a fitted model.

simulate_dynamic_poisson(), simulate_dynamic_binomial()

Simulate data.

Plotting

plot_latent(), plot_fitted(), plot_forecast(), plot_zero_inflation().

Latent dynamics

The latent_dynamics argument of fit_dynamic_model() selects the GMRF state evolution z_t = \mu + \rho z_{t-1} + \varepsilon_t:

"rw"

A first-order random walk, i.e. \rho = 1 fixed (the default; \rho is not sampled).

"ar1"

A stationary AR(1) process, always with an intercept (include_mu = TRUE).

Both share the precision P = D_\rho^\top \mathrm{diag}(1/\sigma^2) D_\rho, where D_\rho is the generalised first-difference operator; it is symmetric tridiagonal, so its bands are formed and the latent full conditionals evaluated in O(N). A scalar \mu (set include_mu = TRUE) adds a drift (RW) or intercept (AR(1)) to the state equation; it enters as a linear term in the latent full conditionals, leaving the sparse precision unchanged.

The otherwise-improper GMRF is anchored by a proper, fixed N(\code{init\_mean}, \code{init\_var}) prior on the first latent state (default N(0, 100)), under both dynamics. Because this prior does not involve (\mu, \rho), the joint full conditional of the drift/intercept and the AR(1) coefficient is conjugate Gaussian, and both are sampled by exact Gibbs draws – \rho from its marginal truncated to the stationary region (-1, 1), then \mu given \rho – with no Metropolis-Hastings correction.

The four innovation structures

The innovations argument controls the distribution of the latent increments \varepsilon_t = z_t - \mu - \rho z_{t-1} (with \rho = 1 for the random walk and \mu = 0 unless a drift/intercept is included):

"gaussian"

\varepsilon_t \sim N(0, \sigma^2) with a single, constant variance.

"t"

A Student-t scale mixture: \varepsilon_t \sim t_\nu(0, \sigma^2), robust to occasional large jumps. The degrees of freedom \nu are sampled.

"mixture"

A finite scale mixture of normals with H components, a flexible heavy-tailed / heteroskedastic increment distribution.

"sv"

Stochastic volatility: \log\sigma^2_t follows an AR(1) process (delegated to the stochvol package).

Zero inflation

Set zeros = "inflated" (or zero_inflation = TRUE) to fit time-constant zero inflation with either family: the observed count is y_t = v_t \tilde y_t, where the gate v_t \sim \mathrm{Bernoulli}(\pi_{\mathrm{open}}) decides whether an observed zero is structural (gate closed) or a sampling zero produced by the Poisson/binomial process (gate open). The gate-open probability \pi_{\mathrm{open}} is a single parameter that does not vary over time or with covariates. See structural_zero_prob().

Every fit stores two flavours of in-sample fitted values and replicates: unconditional draws (fitted, yrep), which include the gate and are the quantities to compare with observed data (posterior predictive checks), and conditional-on-gate-open draws (fitted_open, yrep_open), which come straight from the observation model and describe the latent intensity process. Without zero inflation the pairs coincide exactly. Response forecasts (forecast_y) are always unconditional – the gate is applied. See predict.dynamic_fit().

Forecasting

Forecasts are produced during sampling: fit with horizon = H and H latent states are appended to the state vector and sampled as missing values inside the MCMC (no likelihood contribution, drawn from their GMRF full conditionals). A response is drawn from the observation model at each sampled forecast state. Retrieve the stored draws with forecast.dynamic_fit().

Author(s)

Maintainer: Gregor Zens zens@iiasa.ac.at

References

Zens, G. and Bijak, J. (2026). Dynamic Count Models with Flexible Innovation Processes for Irregular Maritime Migration. The Annals of Applied Statistics. doi:10.1214/26-AOAS2171.

Specify priors for a dynamic count / binomial model

Description

Builds the prior hyperparameters used by fit_dynamic_model(). Called with no arguments it returns weakly informative defaults.

Usage

dynamic_prior(
  var_shape = 2.5,
  var_rate = 0.5,
  df_min = 3,
  df_mean_excess = 6,
  mix_components = 2,
  mix_concentration = 1,
  sv_prior = NULL,
  zi_open_a = 1,
  zi_open_b = 1,
  ar_rho_mean = 0,
  ar_rho_sd = 1,
  mu_mean = 0,
  mu_sd = 1,
  init_mean = 0,
  init_var = 100
)

Arguments

Details

The model places a GMRF latent process on the series z_t (a log-rate for the Poisson family, a logit for the binomial family): either a first-order random walk (latent_dynamics = "rw") or an AR(1) process (latent_dynamics = "ar1"). The increments \varepsilon_t = z_t - \mu - \rho\, z_{t-1} (with \rho = 1 for the random walk and \mu = 0 unless a drift/intercept is included) are given one of four innovation distributions, each governed by some of the priors below.

Innovation variance (all structures). The baseline increment variance \sigma^2 (or, for the "t"/"mixture" structures, the overall scale) has an inverse-gamma prior

\sigma^2 \sim \mathrm{InvGamma}(\code{var\_shape}, \code{var\_rate}).

Smaller var_rate gives a smoother (more strongly shrunk) latent path.

Student-t degrees of freedom (innovations = "t"). The degrees of freedom are modelled as \nu = \code{df\_min} + E, where E \sim \mathrm{Exp}(\mathrm{rate} = 1/\code{df\_mean\_excess}). Thus \nu \ge \code{df\_min} and its prior mean is \code{df\_min} + \code{df\_mean\_excess}. Large \nu approaches the Gaussian case.

Scale mixture (innovations = "mixture"). The mixture uses mix_components variance components, each with an \mathrm{InvGamma}(\code{var\_shape}, \code{var\_rate}) prior, and symmetric Dirichlet weights with concentration mix_concentration.

Stochastic volatility (innovations = "sv"). The log-variance AR(1) priors are delegated to stochvol. Pass a prior specification created with stochvol::specify_priors() via sv_prior, or leave it NULL to use the stochvol defaults.

Zero inflation (zeros = "inflated"). The probability that the observation "gate" is open (i.e. that a zero is an ordinary sampling zero rather than a structural zero) has a \mathrm{Beta}(\code{zi\_open\_a}, \code{zi\_open\_b}) prior. The default Beta(1, 1) is uniform.

AR(1) coefficient (latent_dynamics = "ar1"). When the latent state follows z_t = \mu + \rho\, z_{t-1} + \varepsilon_t, the coefficient \rho is given a Gaussian prior \rho \sim \mathrm{N}(\code{ar\_rho\_mean}, \code{ar\_rho\_sd}^2), truncated to the stationary region \rho \in (-1, 1), and is sampled jointly with \mu by an exact conjugate Gibbs draw (see DynCount-package). AR(1) always carries an intercept (include_mu = TRUE). Under latent_dynamics = "rw" the coefficient is fixed at \rho = 1 and this prior is unused.

Drift / intercept (include_mu = TRUE). A scalar \mu in the state equation z_t = \mu + \rho\, z_{t-1} + \varepsilon_t: a drift under the random walk (\rho = 1) and an intercept under AR(1) (where it is always enabled). It is given a Gaussian prior \mu \sim \mathrm{N}(\code{mu\_mean}, \code{mu\_sd}^2). With include_mu = FALSE (random walk only), \mu = 0 and is not sampled.

Initial state. A proper, fixed \mathrm{N}(\code{init\_mean}, \code{init\_var}) prior (default N(0, 100)) anchors the otherwise-improper GMRF on the first latent state, under both "rw" and "ar1" dynamics; keeping it free of (\mu, \rho) is what makes their update conjugate (see DynCount-package). With the diffuse default the initial state is effectively determined by the data.

Value

An object of class "dynamic_prior": a named list of hyperparameters.

See Also

fit_dynamic_model()

Examples

# Defaults
dynamic_prior()

# Smoother latent path and heavier-tailed t innovations
dynamic_prior(var_rate = 0.1, df_min = 2, df_mean_excess = 3)

Fit a Bayesian dynamic count / binomial time-series model

Description

Fits a GMRF state-space model in which a latent trajectory z_t evolves as either a first-order random walk (latent_dynamics = "rw", the default) or a stationary AR(1) process (latent_dynamics = "ar1"), and the observations are linked to it through a Poisson (log link) or binomial (logit link) observation model. See DynCount-package for an overview.

Usage

fit_dynamic_model(
  y,
  family = c("poisson", "binomial"),
  trials = NULL,
  innovations = c("gaussian", "t", "mixture", "sv"),
  latent_dynamics = c("rw", "ar1"),
  include_mu = FALSE,
  zeros = c("none", "inflated", "missing"),
  zero_inflation = FALSE,
  prior = dynamic_prior(),
  nsave = 4000,
  nburn = 1000,
  thin = 1,
  horizon = 0L,
  forecast_trials = NULL,
  forecast_offset = NULL,
  offset = NULL,
  verbose = FALSE,
  seed = NULL
)

Arguments

Details

Zero handling. Under zeros = "inflated" a latent gate decides, for each observed zero, whether it is structural (gate closed) or an ordinary sampling zero produced by the Poisson/binomial process (gate open); the gate-open probability is a single time-constant parameter. See structural_zero_prob().

Forecasting. Forecasts are produced during sampling: set horizon = H (a positive integer) and H extra latent states are appended to the state vector and sampled as missing values inside the MCMC. A response is drawn from the observation model at each sampled forecast state. The stored forecast draws are retrieved with forecast.dynamic_fit(). With the default horizon = 0 no forecast is produced.

Value

An object of class "dynamic_fit": a list with three elements, each holding a distinct kind of content.

draws

Posterior draws only – a list of matrices/vectors with one row (or element) per stored draw:

z, sig2

Latent states and increment variances.

fitted, yrep

Unconditional fitted means and posterior predictive replicates; under zeros = "inflated" they include the zero-inflation gate.

fitted_open, yrep_open

Their conditional-on-gate-open counterparts: the latent-implied mean, and a replicate drawn straight from the observation model.

gate, pi_open

Zero-inflation gate indicators and gate-open probability. For fits without zero inflation these are stored as placeholders: gate is the constant indicator y > 0 and pi_open is 1.

rho, mu, nu, innov_var

AR(1) coefficient, drift/intercept, Student-t degrees of freedom, and a representative innovation variance whose definition depends on innovations: the estimated constant increment variance \sigma^2 for "gaussian"; the marginal increment variance for "t" (\sigma^2 \nu / (\nu - 2)) and "mixture" (\sigma^2 \sum_h \eta_h \sigma_h^2); and the average of the per-increment variances \exp(h_t) over the series for "sv". For "t", if a draw has \nu \le 2 (possible only when df_min <= 2), the undefined variance factor is replaced by the convention 3.

forecast_z, forecast_y

Latent and response forecasts when horizon >= 1 (NULL otherwise); forecast_y is unconditional, i.e. includes the gate.

Use yrep, not yrep_open, for posterior predictive checks; without zero inflation the conditional and unconditional pairs are identical. See predict.dynamic_fit().

data

The observed inputs: y, trials, the series length n, and the resolved per-observation offset.

spec

The model and MCMC specification: family, innovations, latent_dynamics, include_mu, zeros, prior, nsave, nburn, thin, horizon, and the fixed forecast-period inputs forecast_offset / forecast_trials (NULL when not applicable).

See Also

dynamic_prior(), forecast.dynamic_fit(), plot_fitted(), structural_zero_prob()

Examples

set.seed(1)
sim <- simulate_dynamic_poisson(n = 60, sigma = 0.2, log_rate0 = 2)
fit <- fit_dynamic_model(sim$y, family = "poisson", nsave = 300, nburn = 200)
summary(fit)

Generic forecast function

Description

Generic for extracting forecasts from a fitted model. See forecast.dynamic_fit() for the method for "dynamic_fit" objects.

Usage

forecast(object, ...)

Arguments

Value

The value returned by the method dispatched on object; for "dynamic_fit" objects a "dynamic_forecast" object (see forecast.dynamic_fit()).

Forecast a fitted dynamic model

Description

Returns the forecast draws that were produced during MCMC sampling. Forecasts are generated inside the sampler when the model is fitted with horizon = H >= 1. This method extracts and summarises those stored draws.

Usage

## S3 method for class 'dynamic_fit'
forecast(object, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Details

Zero inflation. For fits with zeros = "inflated" the response forecast is unconditional: each draw includes the zero-inflation gate, so structural zeros are reproduced and the draws are directly comparable to future observations (like the in-sample yrep, unlike yrep_open). The latent forecast summary in ⁠$latent⁠ is gate-free by construction.

Value

An object of class "dynamic_forecast": a list with summary (one row per horizon 1, \dots, H, on the response scale), latent (summary of the forecast latent path), draws (the predictive draws, draws x H), final (the single-row summary of the final H-step-ahead forecast), final_draws (the predictive draws at horizon H), and the forecast horizon horizon.

Examples

sim <- simulate_dynamic_poisson(60, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 300, nburn = 200, horizon = 8)
fc <- forecast(fit)
fc$summary

Weekly Mediterranean crossings (Mediterranean example)

Description

A longer weekly count series of irregular sea crossings on the Mediterranean route, beginning in autumn 2015. The counts are larger in magnitude with fewer zeros than uk_weekly, which makes the series useful for the plain Poisson model and the robust (Student-t) and stochastic-volatility innovation structures. The series is used in the irregular-migration application of Zens and Bijak (2026).

Usage

med_weekly

Format

A data frame with 494 rows and 2 variables:

week

ISO week label, e.g. "2015-W40".

count

Non-negative integer count of weekly crossings.

Source

Weekly aggregates of detected irregular Mediterranean crossings, compiled from operational/agency records as described in Zens and Bijak (2026), doi:10.1214/26-AOAS2171.

References

Zens, G. and Bijak, J. (2026). Dynamic Count Models with Flexible Innovation Processes for Irregular Maritime Migration. The Annals of Applied Statistics. doi:10.1214/26-AOAS2171.

Examples

data(med_weekly)
summary(med_weekly$count)

Plot method for fitted dynamic models

Description

A convenience wrapper that dispatches to the dedicated plotting functions.

Usage

## S3 method for class 'dynamic_fit'
plot(x, which = c("fitted", "latent", "forecast", "zeros"), ...)

Arguments

Value

Invisibly, the result of the underlying plotting function.

Plot observed versus fitted values

Description

Observed series with the posterior median fitted mean and a credible band on the response scale.

Usage

plot_fitted(object, level = 0.95, ...)

Arguments

Value

Invisibly, the fitted summary.

Examples

sim <- simulate_dynamic_poisson(60, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 300, nburn = 200)
plot_fitted(fit)

Plot observed history, fitted values and forecast with uncertainty

Description

Shows the observed series together with the in-sample fitted mean (median and credible band) and the forecast (median and credible band) over the forecast horizon. The forecast draws are those produced during sampling, so the model must have been fitted with horizon >= 1 (see forecast.dynamic_fit()).

Usage

plot_forecast(object, level = 0.95, ...)

Arguments

Value

Invisibly, the "dynamic_forecast" object.

Examples

sim <- simulate_dynamic_poisson(60, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 300, nburn = 200, horizon = 12)
plot_forecast(fit)

Plot the fitted latent trajectory

Description

Shows the posterior median of the latent state (log-rate for the Poisson family, logit for the binomial family) with a credible band.

Usage

plot_latent(object, level = 0.95, ...)

Arguments

Value

Invisibly, the summary data frame used for plotting.

Examples

sim <- simulate_dynamic_poisson(60, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 300, nburn = 200)
plot_latent(fit)

Plot zero-inflation diagnostics

Description

Bar chart of the posterior probability that each observed zero is a structural zero. Requires a model fitted with zeros = "inflated".

Usage

plot_zero_inflation(object, ...)

Arguments

Value

Invisibly, the structural_zero_prob() data frame.

Examples

sim <- simulate_dynamic_poisson(80, 0.2, 2, zero_inflation = 0.3, seed = 1)
fit <- fit_dynamic_model(sim$y, zero_inflation = TRUE, nsave = 300, nburn = 200)
plot_zero_inflation(fit)

In-sample fitted values and posterior predictive replicates

Description

In-sample fitted values and posterior predictive replicates

Usage

## S3 method for class 'dynamic_fit'
predict(
  object,
  type = c("mean", "response"),
  conditional = FALSE,
  probs = c(0.025, 0.5, 0.975),
  ...
)

Arguments

Value

A list with summary (a data frame, one row per observation) and draws (the underlying draws matrix, draws x time).

Examples

sim <- simulate_dynamic_poisson(40, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 200, nburn = 100)
head(predict(fit)$summary)

Simulate a binomial dynamic series

Description

Generates a latent logit process (random walk or AR(1)) and binomial counts, optionally with structural (zero-inflation) zeros.

Usage

simulate_dynamic_binomial(
  n,
  sigma,
  trials,
  logit0 = 0,
  zero_inflation = 0,
  rho = 1,
  mu = 0,
  offset = 0,
  seed = NULL
)

Arguments

Value

A list with components y (successes), trials, logit (latent path z_t), prob (plogis(offset + logit)), offset, and structural (logical, TRUE for structural zeros).

Examples

sim <- simulate_dynamic_binomial(n = 50, sigma = 0.15, trials = 40, seed = 1)
head(sim$y)
# with structural zeros:
zi <- simulate_dynamic_binomial(50, 0.15, trials = 40, zero_inflation = 0.2, seed = 1)
mean(zi$structural)

Simulate a Poisson dynamic series

Description

Generates a latent log-rate process (random walk or AR(1)) and Poisson counts, optionally with zero inflation.

Usage

simulate_dynamic_poisson(
  n,
  sigma,
  log_rate0 = 1,
  zero_inflation = 0,
  rho = 1,
  mu = 0,
  offset = 0,
  seed = NULL
)

Arguments

Value

A list with components y (observed counts), log_rate (the latent log-rate path z_t), rate (the mean exp(offset + log_rate)), offset, and structural (logical, TRUE where a structural zero was forced).

Examples

sim <- simulate_dynamic_poisson(n = 50, sigma = 0.2, log_rate0 = 2,
                           zero_inflation = 0.2, seed = 1)
table(sim$y == 0, sim$structural)

Posterior probability that each observed zero is structural

Description

For a zero-inflated fit (Poisson or binomial), returns for every observed zero the posterior probability that it is a structural zero (gate closed) rather than an ordinary sampling zero generated by the observation process. The zero inflation is time-constant: a single gate-open probability governs all observations (it does not vary over time or with covariates).

Usage

structural_zero_prob(object, zeros_only = TRUE)

Arguments

Details

The model introduces a latent gate indicator v_t \in \{0, 1\}: when the gate is open (v_t = 1) the count comes from the observation model (Poisson or binomial); when closed (v_t = 0) the count is a structural zero. The sampler stores v_t for each draw, so the structural-zero probability is

\Pr(\text{structural} \mid y_t = 0) = 1 - \overline{v_t},

the posterior mean of the gate being closed. Non-zero observations are always sampling observations and have structural probability 0.

Value

A data frame with columns time, observed, p_structural (posterior probability the zero is structural) and p_sampling (1 - p_structural).

Examples

sim <- simulate_dynamic_poisson(60, 0.2, 2, zero_inflation = 0.25, seed = 1)
fit <- fit_dynamic_model(sim$y, zero_inflation = TRUE, nsave = 300, nburn = 200)
structural_zero_prob(fit)

Summarise a fitted dynamic model

Description

Summarise a fitted dynamic model

Usage

## S3 method for class 'dynamic_fit'
summary(object, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Value

An object of class "summary.dynamic_fit" containing posterior summaries of the global parameters (the innovation standard deviation innov_sd – the estimated SD for "gaussian" innovations, the marginal SD for "t" and "mixture", and the series-average SD for "sv" – and where relevant the t degrees of freedom, the AR(1) coefficient, the drift/intercept, and the zero-inflation gate-open probability) and the fitted-value summary.

Examples

sim <- simulate_dynamic_poisson(50, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 200, nburn = 100)
summary(fit)

Weekly English Channel crossings (UK example)

Description

A weekly count series of irregular small-boat crossings of the English Channel towards the United Kingdom, beginning in early 2018. The series has frequent zeros in its early weeks, which makes it a useful example for zero inflation with the Poisson model. The series is used in the irregular-migration application of Zens and Bijak (2026).

Usage

uk_weekly

Format

A data frame with 376 rows and 2 variables:

week

ISO week label, e.g. "2018-W01".

count

Non-negative integer count of weekly crossings.

Source

Weekly aggregates of detected irregular English Channel crossings, compiled from operational/agency records as described in Zens and Bijak (2026), doi:10.1214/26-AOAS2171.

References

Zens, G. and Bijak, J. (2026). Dynamic Count Models with Flexible Innovation Processes for Irregular Maritime Migration. The Annals of Applied Statistics. doi:10.1214/26-AOAS2171.

Examples

data(uk_weekly)
plot(uk_weekly$count, type = "h")
mean(uk_weekly$count == 0)