Package {DynCount}

DynCount: Bayesian Dynamic Models for Count Time Series

Description

The DynCount package fits 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 a Poisson (log link), a binomial (logit link) or a multinomial (additive-log-ratio 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, per-category ALR shifts for the multinomial shares). This 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). It supports several innovation structures (Gaussian, Student-t, finite scale mixture, stochastic volatility) and, for the Poisson and binomial families, zero inflation with a time-constant gate-open probability.

Main entry points

fit_dynamic_model()

Fit a Poisson, binomial or multinomial dynamic model (random walk or stationary AR(1)).

forecast() / predict()

Posterior predictive forecasts for any horizon (by forward simulation from the posterior draws), or the in-sample fitted values and replicates.

summary()

Posterior summaries of a fitted model.

simulate_dynamic_poisson(), simulate_dynamic_binomial(), simulate_dynamic_multinomial()

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).

"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.

Multinomial (choice-count) series

With family = "multinomial" the data are an n \times K matrix of category counts with known row totals N_t. One category b (by default the one with the largest total count) is the baseline, and each of the other K - 1 categories has its own latent additive-log-ratio series z_{t,k} = \log(p_{t,k} / p_{t,b}), so that p_t is the softmax of (z_{t,\cdot}, 0) and y_t \sim \mathrm{Multinomial}(N_t, p_t). Every ALR series follows the selected latent dynamics with the selected innovation structure, but the series share no parameters. Each series has its own innovation variance (and auxiliary innovation parameters), its own \rho and \mu, and its own copy of the prior. They are coupled only through the multinomial likelihood, and each ALR series is updated in turn with the other categories held at their current values. With K = 2 and the second column as baseline the model coincides exactly with the binomial model. The model is not invariant to the choice of baseline (the dynamics are placed on the log-ratios relative to it), so a large category with a stable share is the natural choice. Zero inflation is not available for this family currently. Rows with N_t = 0 are uninformative and handled as missing. Posterior draws carry a trailing category dimension, see fit_dynamic_model().

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 estimated from the data.

"mixture"

A finite scale mixture of normals with mix_components components (see dynamic_prior()), resulting in a flexible 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 zero inflation with the Poisson or the binomial 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 are always unconditional, with the gate applied to each forecast draw. See predict.dynamic_fit().

Forecasting

Future latent states carry no likelihood, so their posterior given a draw of the parameters and of the last in-sample state is the transition law of the latent process. forecast.dynamic_fit() therefore forward-simulates the latent path from every stored posterior draw (with increments from the fitted innovation structure) and draws a response from the observation model at each simulated state. Forecasts can be requested for any horizon after fitting, e.g. forecast(fit, horizon = 8). Fitting with horizon = H stores an H-step forecast in the 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, 20(2), 1671–1690. doi:10.1214/26-AOAS2171.

Specify priors for a dynamic 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 = 0.01,
  var_rate = 0.01,
  df_min = 3,
  df_mean_excess = 6,
  mix_components = 2,
  mix_concentration = 1,
  mix_var_shape = 2.5,
  mix_var_rate = 0.5,
  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, and one additive-log-ratio series per non-baseline category for the multinomial 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}).

The default \mathrm{InvGamma}(0.01, 0.01) is weakly informative for typical increment standard deviations, but it is not scale-free. For very smooth or short series, with increment standard deviations of a few hundredths, the results can be sensitive to this prior. In that regime, check the sensitivity of the results to a smaller var_rate (which lets \sigma become smaller) and expect slower mixing, because a nearly constant latent path and a small \sigma are strongly dependent a posteriori. In the multinomial family the same prior is applied independently to every non-baseline category.

Mixture component variances (innovations = "mixture"). The relative variances of the mixture components multiply the overall scale and are given their own \mathrm{InvGamma}(\code{mix\_var\_shape}, \code{mix\_var\_rate}) prior (default \mathrm{InvGamma}(2.5, 0.5)). They are kept moderately informative on purpose, as with a vague prior an empty component would be drawn from an extremely heavy-tailed distribution and the split between the overall scale and the component variances is only weakly identified.

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{mix\_var\_shape}, \code{mix\_var\_rate}) prior (see above), and symmetric Dirichlet weights with concentration mix_concentration.

Stochastic volatility (innovations = "sv"). The log-variance h_t of the increments follows h_t = \mu_h + \phi (h_{t-1} - \mu_h) + \sigma_h \eta_t, and its 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). 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. 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()

# An informative variance prior and heavier-tailed t innovations
dynamic_prior(var_shape = 2.5, var_rate = 0.5, df_min = 2, df_mean_excess = 3)

Fit a Bayesian dynamic count / binomial / multinomial 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), binomial (logit link) or multinomial (additive-log-ratio link) observation model. See DynCount-package for an overview.

Usage

fit_dynamic_model(
  y,
  family = c("poisson", "binomial", "multinomial"),
  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,
  baseline = "largest",
  verbose = FALSE,
  seed = NULL
)

Arguments

Details

Estimation. The model is estimated by Metropolis-within-Gibbs MCMC. The latent states are updated one site at a time by adaptive random-walk Metropolis steps that use their Gaussian Markov random field (GMRF) full conditionals. States without an observation (the initial state, structural zeros, zero-total rows) are drawn exactly from their Gaussian full conditionals. The innovation parameters, the drift/intercept \mu and the AR(1) coefficient \rho are updated by Gibbs steps (with a Metropolis step for the Student-t degrees of freedom).

Multinomial family. With family = "multinomial", y is an n \times K matrix of category counts (rows = time, columns = categories, K \ge 2); the row totals N_t are treated as known trials. One category b is the baseline, and the remaining K - 1 categories each get their own latent additive-log-ratio (ALR) series z_{t,k} = \log(p_{t,k} / p_{t,b}), so that

y_t \sim \mathrm{Multinomial}(N_t, p_t), \qquad p_{t,k} = \frac{\exp(o_{t,k} + z_{t,k})}{1 + \sum_{j \ne b} \exp(o_{t,j} + z_{t,j})}, \qquad p_{t,b} = \frac{1}{1 + \sum_{j \ne b} \exp(o_{t,j} + z_{t,j})},

with known offsets o_{t,k} (zero by default). The chosen latent dynamics, innovation structure and drift/intercept setting apply to every ALR series, but no parameters are shared across categories. Each series has its own innovation variance (and, where relevant, degrees of freedom, mixture components or volatility path), its own \rho and \mu, and its own copy of prior. The series are coupled only through the multinomial likelihood, and each series is updated with the other categories held at their current values. A row with N_t = 0 carries no information about the shares and is handled like a missing observation. Zero inflation is not available for this family (zeros must be "none"). Running time grows linearly in K - 1.

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 parameter that is constant over time. See structural_zero_prob(). Under zeros = "missing" the observed zeros are treated as missing values.

Forecasting. Forecasts are obtained by forward simulation from the posterior draws, so they can be computed after fitting for any horizon with forecast.dynamic_fit(). Setting horizon = H additionally simulates an H-step forecast right after sampling and stores it in the fit (under the same seed).

Reproducibility. With a seed, the sampler and the fit-time forecast run with that seed, and the previous state of the global random number generator is restored afterwards.

Value

An object of class "dynamic_fit": a list with three elements.

draws

Posterior draws only. For the Poisson and binomial families a list of matrices/vectors with one row (or element) per stored draw:

z

Latent states aligned with the observations (⁠draws x n⁠).

z0

The initial latent state, one period before the first observation, which carries the init_mean / init_var prior.

sig2

Variance of the increment leading into each z_t (⁠draws x n⁠); the first column is the increment from z0.

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; NULL unless zeros = "inflated".

rho, mu

AR(1) coefficient (1 under the random walk) and drift/intercept (0 unless include_mu = TRUE).

innov_var

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 in-sample transitions 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.

scale

The overall innovation scale \sigma^2 (NULL for "sv").

nu

Student-t degrees of freedom (NULL unless innovations = "t").

mix_weight, mix_var

Mixture weights and component variances (\sigma^2 \sigma_h^2), ⁠draws x mix_components⁠ (NULL unless innovations = "mixture"). The components are exchangeable and not identified individually. Within each draw they are stored in order of increasing variance, which is a labelling convention rather than an identification, so per-component summaries are meaningful only if the components are clearly separated. Label-invariant quantities, such as innov_var and the forecasts, are unaffected.

sv_mu, sv_phi, sv_sigma

Level, persistence and volatility of the AR(1) log-variance process (NULL unless innovations = "sv").

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().

For the multinomial family the same components carry an extra trailing category dimension, named by the category labels: latent quantities (z, sig2, forecast_z, mix_weight, mix_var) are arrays with one slice per non-baseline category, and z0, rho, mu, innov_var, scale, nu and the ⁠sv_*⁠ parameters are ⁠draws x (K - 1)⁠ matrices; response-scale quantities (fitted – the expected counts N_t p_{t,k} –, yrep, forecast_y, and the additional fitted_prob / forecast_prob holding the category shares) are ⁠draws x time x K⁠ arrays including the baseline, in the original column order. fitted_open, yrep_open, gate and pi_open are NULL.

data

The observed inputs: y, trials (the row totals for the multinomial family), the series length n, and the resolved per-observation offset. For the multinomial family also K, the categories (column labels), and the baseline index and baseline_name.

spec

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

See Also

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

Examples

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

# multinomial choice counts: three categories, baseline chosen automatically
simm <- simulate_dynamic_multinomial(n = 40, sigma = 0.15, trials = 100,
                                     alr0 = c(-1, -0.5), seed = 2)
fitm <- fit_dynamic_model(simm$y, family = "multinomial",
                          nsave = 200, nburn = 100, seed = 2)
fitm
summary(fitm)$params

Forecast a fitted dynamic model

Description

Returns posterior predictive forecasts for the H periods after the end of the series. The forecasts are obtained by forward simulation from the stored posterior draws, so any horizon can be requested after fitting. If the model was fitted with horizon >= 1, the forecast draws stored in the fit are returned unless horizon, forecast_offset or forecast_trials is supplied.

Usage

## S3 method for class 'dynamic_fit'
forecast(
  object,
  horizon = NULL,
  forecast_offset = NULL,
  forecast_trials = NULL,
  probs = c(0.025, 0.5, 0.975),
  seed = NULL,
  ...
)

Arguments

Details

For every stored posterior draw, the latent path is propagated from the last in-sample state z_n with z_{n+h} = \mu + \rho z_{n+h-1} + \varepsilon_{n+h}, drawing the increments from the fitted innovation structure (Gaussian with the drawn variance, Student-t with the drawn scale and degrees of freedom, the drawn scale mixture, or the stochastic-volatility process continued from its last in-sample log-variance). A response is then drawn from the observation model at each simulated state. The intervals therefore reflect parameter, state and innovation uncertainty.

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.

Multinomial family. Response forecasts are multinomial draws with the totals given by forecast_trials. The summaries are in long format with a category column. summary and final cover all K categories on the count scale, prob gives the forecast category shares, and latent the K - 1 ALR series. draws is a ⁠draws x H x K⁠ array and final_draws a ⁠draws x K⁠ matrix.

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), latent_draws (the latent forecast draws), final (the single-row summary of the final H-step-ahead forecast), final_draws (the predictive draws at horizon H), the forecast horizon, and the forecast_offset / forecast_trials that were used. For the multinomial family additionally prob (share summary) and prob_draws (⁠draws x H x K⁠ share draws); see Details for the layout.

See Also

fit_dynamic_model(), plot_forecast()

Examples

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

Weekly Mediterranean crossings (Mediterranean example)

Description

A longer weekly count series of irregular sea crossings on the Mediterranean route, covering the ISO weeks 2015-W40 to 2025-W11 (494 weeks). The counts are larger in magnitude with fewer zeros than uk_weekly. 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 3 variables:

week

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

count

Non-negative integer count of weekly crossings.

date

The Monday of the ISO week, as a Date.

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, 20(2), 1671–1690. doi:10.1214/26-AOAS2171.

Examples

summary(med_weekly$count)
plot(med_weekly$date, med_weekly$count, type = "h", xlab = "week", ylab = "crossings")

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. For the multinomial family, one panel per category showing the observed and expected counts N_t p_{t,k}.

Usage

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

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, seed = 1)
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 is the one stored in the fit, or a new one computed with forecast.dynamic_fit() when horizon (or another forecast input) is supplied. For the multinomial family, one panel per category.

Usage

plot_forecast(
  object,
  level = 0.95,
  category = NULL,
  horizon = NULL,
  forecast_offset = NULL,
  forecast_trials = NULL,
  seed = NULL,
  ...
)

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, seed = 1)
plot_forecast(fit, horizon = 12)

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, one additive-log-ratio series per non-baseline category for the multinomial family) with a credible band.

Usage

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

Arguments

Value

Invisibly, the summary data frame used for plotting (in long format with a category column for the multinomial family).

Examples

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

Plot zero-inflation diagnostics

Description

Bar chart of the posterior probability that each observed zero is a structural zero (see structural_zero_prob()), one bar per observed zero labelled by its time index, with a dotted reference line at 0.5. Requires a model fitted with zeros = "inflated" (Poisson or binomial family).

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,
                         seed = 1)
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", "prob"),
  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). For the multinomial family summary is in long format with one row per observation and category (columns time, category, observed, ...) and draws is a ⁠draws x time x K⁠ array.

Examples

sim <- simulate_dynamic_poisson(40, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 200, nburn = 100, seed = 1)
head(predict(fit)$summary)                     # posterior of the mean of y
head(predict(fit, type = "response")$summary)  # posterior predictive

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

forecast

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 multinomial dynamic series

Description

Generates K - 1 independent latent additive-log-ratio (ALR) processes (random walk or AR(1)), one per non-baseline category, and multinomial choice counts with known totals. The ALR series z_{t,k} = \log(p_{t,k} / p_{t,b}) share no parameters: sigma, rho and mu may each be a single value (recycled) or a vector of length K - 1 giving one value per non-baseline category.

Usage

simulate_dynamic_multinomial(
  n,
  sigma,
  trials,
  alr0 = c(0, 0),
  baseline = length(alr0) + 1L,
  rho = 1,
  mu = 0,
  offset = 0,
  categories = NULL,
  seed = NULL
)

Arguments

Value

A list with components y (an ⁠n x K⁠ matrix of counts with the baseline in column baseline), trials, alr (the ⁠n x (K - 1)⁠ latent ALR paths z_{t,k}), prob (the ⁠n x K⁠ matrix of category probabilities), offset, baseline (the column index) and categories.

Examples

sim <- simulate_dynamic_multinomial(n = 50, sigma = 0.15, trials = 200,
                                    alr0 = c(-0.5, 0.5), seed = 1)
head(sim$y)
colSums(sim$y)
# fit on the simulated ALR scale by passing the simulated baseline
fit <- fit_dynamic_model(sim$y, family = "multinomial", baseline = sim$baseline,
                         nsave = 200, nburn = 100, seed = 1)
summary(fit)$params

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. A single gate-open probability governs all observations (it does not vary over time or with covariates), while the gate itself is drawn separately for every observation.

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,
                         seed = 1)
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 (params, one row per parameter) and of the fitted values (fitted). The rows of params and the draws they summarise are

innov_sd

sqrt(draws$innov_var): the estimated increment SD for "gaussian" innovations, the marginal SD for "t" and "mixture", and the root of the series-average variance for "sv".

t_df

draws$nu, the Student-t degrees of freedom.

sv_mu, sv_phi, sv_sigma

draws$sv_mu, draws$sv_phi and draws$sv_sigma, the parameters of the AR(1) log-variance process.

ar1_rho

draws$rho, the AR(1) coefficient.

drift_mu / intercept_mu

draws$mu, the random-walk drift or the AR(1) intercept.

gate_open_prob

draws$pi_open, the probability that the zero-inflation gate is open (one minus the structural-zero probability).

Rows are included only where relevant. For "mixture" innovations only innov_sd is reported. The mixture components are exchangeable and not identified individually, so per-component summaries are not given. The draws of the weights and component variances are in draws$mix_weight and draws$mix_var. For the multinomial family every latent parameter is reported once per non-baseline category, with the category label in square brackets (e.g. innov_sd[B]), and the fitted summary is a long data frame with a category column covering all K categories.

Examples

sim <- simulate_dynamic_poisson(50, 0.2, 2, seed = 1)
fit <- fit_dynamic_model(sim$y, nsave = 200, nburn = 100, seed = 1)
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, covering the ISO weeks 2018-W01 to 2025-W11 (376 weeks). The series has frequent zeros in its early weeks (43% of the first 130 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 3 variables:

week

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

count

Non-negative integer count of weekly crossings.

date

The Monday of the ISO week, as a Date.

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, 20(2), 1671–1690. doi:10.1214/26-AOAS2171.

Examples

plot(uk_weekly$date, uk_weekly$count, type = "h", xlab = "week", ylab = "crossings")
mean(uk_weekly$count == 0)