Dynamic Models for Poisson and Binomial Time Series

Introduction

DynCount fits state-space models to non-Gaussian time series. A latent trajectory $z_t$ evolves with one of two dynamics, \[ z_t = \mu + \rho\, z_{t-1} + \varepsilon_t, \] a first-order random walk (latent_dynamics = "rw", i.e. $\rho = 1$) or a stationary AR(1) process (latent_dynamics = "ar1", with $\rho$ estimated and constrained to $(-1, 1)$). The scalar $\mu$ is zero unless a drift/intercept is included. The observations are linked to it through one of two observation models:

Poisson, with a log link: $y_t \sim \mathrm{Poisson}(e^{z_t})$;
Binomial, with a logit link: $y_t \sim \mathrm{Binomial}(m_t,\, \mathrm{logit}^{-1}(z_t))$, where the trials $m_t$ are known.

An optional known offset $o_t$ may be added to the linear predictor of either observation model – a log-exposure for the Poisson mean, $e^{o_t + z_t}$, or a shift of the binomial logit. It is a fixed, user-supplied input, not part of the latent process $z_t$, and defaults to zero.

The model is estimated by MCMC. The distribution of the increments $\varepsilon_t = z_t - \mu - \rho z_{t-1}$ is controlled by the innovations argument and can be Gaussian, Student-t, a finite scale mixture of normals, or a stochastic volatility process. For either family, zeros can optionally be handled by time-constant zero inflation.

The package implements and extends the methodology of Zens and Bijak (2026), The Annals of Applied Statistics, doi:10.1214/26-AOAS2171.

Note: the MCMC runs below use short chains (nsave = 1000, nburn = 1000) and, for the two shipped series, a shortened window, so that the vignette builds quickly. For real analyses use longer chains and the full series.

Simulating data

The simulation helpers generate data with a known latent path, which is useful for checking recovery.

sim <- simulate_dynamic_poisson(n = 80, sigma = 0.18, log_rate0 = 2.5, seed = 1)
str(sim, max.level = 1)
#> List of 5
#>  $ y         : num [1:80] 10 15 13 9 7 12 15 11 19 18 ...
#>  $ log_rate  : num [1:80] 2.5 2.39 2.42 2.27 2.56 ...
#>  $ rate      : num [1:80] 12.18 10.88 11.25 9.68 12.9 ...
#>  $ offset    : num [1:80] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ structural: logi [1:80] FALSE FALSE FALSE FALSE FALSE FALSE ...

plot(sim$y, type = "h", xlab = "time", ylab = "count",
     main = "Simulated Poisson random walk")
lines(sim$rate, col = "steelblue", lwd = 2)

Fitting a Poisson model (no zero inflation)

The main entry point is fit_dynamic_model(). The defaults give an ordinary Poisson random walk with Gaussian increments.

fit <- fit_dynamic_model(sim$y, family = "poisson",
                         nsave = NSAVE, nburn = NBURN, seed = 1)
fit
#> <dynamic_fit>
#>   family      : poisson
#>   dynamics    : rw  (rho = 1)
#>   innovations : gaussian
#>   zeros       : none
#>   observations: 80  (zeros: 0)
#>   draws kept  : 1000

summary(fit)
#> Dynamic count model summary
#>   family = poisson | dynamics = rw | innovations = gaussian | zeros = none
#>   80 observations (0 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>            mean     sd   q2.5    q50  q97.5
#> innov_sd 0.2207 0.0264 0.1777 0.2189 0.2798
#> 
#> Fitted values: range of posterior means [10.56, 97.21]

plot_fitted() overlays the posterior of the fitted mean on the data, and plot_latent() shows the latent log-rate trajectory with a credible band.

plot_fitted(fit)

plot_latent(fit)

Forecasting

Forecasts are produced during sampling. Request the horizon when fitting with horizon = H: H extra latent states are appended to the state vector and sampled as missing values inside the MCMC, and observations are drawn from the observation model at the sampled forecast states, so the intervals propagate parameter, state and innovation uncertainty. forecast() then extracts and summarises the stored draws.

fit_fc <- fit_dynamic_model(sim$y, family = "poisson", horizon = 8,
                            nsave = NSAVE, nburn = NBURN, seed = 1)
fc <- forecast(fit_fc)
fc                                 # prints the forecast path
#> <dynamic_forecast> poisson, horizon 8
#>  horizon   mean       sd q2.5  q50   q97.5
#>        1 61.550 17.99331   34 59.0 105.000
#>        2 63.407 23.32907   29 60.0 121.000
#>        3 64.957 28.61248   27 60.0 136.025
#>        4 65.246 33.52876   22 56.0 148.050
#>        5 67.330 39.70758   21 57.0 170.050
#>        6 69.450 45.93423   17 56.0 181.000
#>        7 71.258 51.30288   16 56.5 209.025
#>        8 72.124 55.56948   15 56.0 220.050

fc$final                           # the single 8-step-ahead forecast
#>   horizon   mean       sd q2.5 q50  q97.5
#> 8       8 72.124 55.56948   15  56 220.05

plot_forecast(fit_fc)

The object stores the full forecast path (fc$summary, one row per horizon) and, separately, the final h-step-ahead prediction (fc$final, fc$final_draws). For long horizons, or with the heavier-tailed innovation structures ("t", "mixture", "sv"), use longer chains than the short ones shown here and, as with any MCMC output, check convergence before relying on the results.

AR(1) latent dynamics

Setting latent_dynamics = "ar1" estimates an autoregressive coefficient $\rho$ instead of fixing it at 1, jointly with an intercept $\mu$. The pair is drawn by an exact conjugate Gibbs step, with $\rho$ truncated to the stationary region $(-1, 1)$, so the posterior places mass only on stationary processes. The random walk is the special case $\rho = 1$. AR(1) always includes an intercept: include_mu is enabled automatically, giving the process a non-zero stationary mean $\mu / (1 - \rho)$.

# a genuinely stationary AR(1) log-rate: stationary mean 4, so mu = 4 * (1 - rho)
sim_ar <- simulate_dynamic_poisson(150, sigma = 0.2, log_rate0 = 4,
                                   rho = 0.9, mu = 0.4, seed = 3)
# no need to set include_mu: AR(1) enables the intercept automatically
fit_ar <- fit_dynamic_model(sim_ar$y, latent_dynamics = "ar1",
                            nsave = NSAVE, nburn = NBURN, seed = 3)
summary(fit_ar)                    # reports the posteriors of ar1_rho and mu
#> Dynamic count model summary
#>   family = poisson | dynamics = ar1 | innovations = gaussian | zeros = none
#>   150 observations (0 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>                mean     sd   q2.5    q50  q97.5
#> innov_sd     0.2069 0.0169 0.1777 0.2057 0.2433
#> ar1_rho      0.8802 0.0456 0.7858 0.8799 0.9634
#> intercept_mu 0.4725 0.1819 0.1340 0.4685 0.8482
#> 
#> Fitted values: range of posterior means [20.80, 121.18]

Offset

For Poisson data with varying exposure, pass a known offset (a log-exposure term); the mean becomes $\exp(\text{offset}_t + z_t)$. Supply forecast_offset at fitting time for the future exposures.

expo  <- log(runif(120, 50, 200))  # known exposure, e.g. population at risk
sim_o <- simulate_dynamic_poisson(120, sigma = 0.12, log_rate0 = -3.5,
                                  offset = expo, seed = 4)
fit_o <- fit_dynamic_model(sim_o$y, offset = expo, horizon = 6,
                           forecast_offset = log(120),
                           nsave = NSAVE, nburn = NBURN, seed = 4)
forecast(fit_o)$final
#>   horizon   mean       sd q2.5 q50  q97.5
#> 6       6 12.906 8.635725    2  11 34.025

A robust (Student-t) fit on the example data

The package ships two real weekly count series of irregular maritime crossings. med_weekly (Mediterranean route) has larger counts and few zeros; the Student-t innovation makes the latent path robust to the occasional large week-to-week jump. For a fast build we demonstrate this using a recent window of the series.

data(med_weekly)
med <- tail(med_weekly$count, 120)
fit_med <- fit_dynamic_model(med, family = "poisson",
                             innovations = "t",
                             nsave = NSAVE, nburn = NBURN, seed = 2)
summary(fit_med)
#> Dynamic count model summary
#>   family = poisson | dynamics = rw | innovations = t | zeros = none
#>   120 observations (1 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>            mean     sd   q2.5    q50  q97.5
#> innov_sd 1.6216 0.1671 1.3348 1.6080 1.9943
#> t_df     4.0197 1.0700 3.0062 3.6499 7.3360
#> 
#> Fitted values: range of posterior means [1.39, 14133.51]

The posterior of the degrees-of-freedom parameter t_df indicates how heavy the increment tails are, with smaller values indicating heavier tails.

Zero inflation and structural zeros

uk_weekly (English Channel crossings) has many zeros. Turning on zero inflation lets the model separate structural zeros (a time-constant gate that switches the count off) from sampling zeros (zeros the Poisson process produces on its own). We use the zero-heavy early window of the series here.

data(uk_weekly)
uk <- uk_weekly$count[1:130]
mean(uk == 0)                         # many zeros
#> [1] 0.4307692

fit_zip <- fit_dynamic_model(uk, family = "poisson",
                             zero_inflation = TRUE,
                             nsave = NSAVE, nburn = NBURN, seed = 3)
summary(fit_zip)
#> Dynamic count model summary
#>   family = poisson | dynamics = rw | innovations = gaussian | zeros = inflated
#>   130 observations (56 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>                  mean     sd   q2.5    q50  q97.5
#> innov_sd       0.8658 0.0848 0.7274 0.8583 1.0552
#> gate_open_prob 0.6974 0.0484 0.6012 0.6993 0.7909
#> 
#> Fitted values: range of posterior means [0.11, 190.07]

structural_zero_prob() reports, for each observed zero, the posterior probability that it is structural.

sz <- structural_zero_prob(fit_zip, zeros_only = TRUE)
head(sz, 10)
#>    time observed p_structural p_sampling
#> 1     1        0        0.483      0.517
#> 2     2        0        0.498      0.502
#> 3     3        0        0.562      0.438
#> 4     4        0        0.698      0.302
#> 5     6        0        0.629      0.371
#> 6     7        0        0.472      0.528
#> 7     8        0        0.412      0.588
#> 8     9        0        0.391      0.609
#> 9    10        0        0.367      0.633
#> 10   11        0        0.360      0.640

plot_zero_inflation(fit_zip)

Interpreting the table: p_structural close to 1 flags a zero that the latent rate cannot easily explain (e.g., the underlying rate was high, so a Poisson zero would be unlikely), whereas p_structural near 0 marks a zero that is consistent with a genuinely low rate.

Conditional versus unconditional fits and replicates

Under zero inflation the observed count is $y_t = v_t \tilde y_t$, where the gate $v_t \sim \mathrm{Bernoulli}(\pi_{\text{open}})$ switches the count off and $\tilde y_t$ comes from the Poisson/binomial observation model. The fit stores both flavours of in-sample quantities:

Unconditional (fitted, yrep; the defaults returned by predict()): the gate is included, so the replicates reproduce the structural zeros. Use these for posterior predictive checks.
Conditional on the gate being open (fitted_open, yrep_open; predict(fit, conditional = TRUE)): the latent-implied mean and a replicate drawn straight from the observation model, describing the latent intensity process.

For models without zero inflation the two versions are identical. Response forecasts (forecast(), forecast_y) are always unconditional – the gate is applied to each forecast draw.

# zero proportion in the data vs both flavours of replicate
c(observed      = mean(uk == 0),
  yrep          = mean(fit_zip$draws$yrep == 0),        # gate applied: comparable
  yrep_open     = mean(fit_zip$draws$yrep_open == 0))   # gate-open only: too few zeros
#>  observed      yrep yrep_open 
#> 0.4307692 0.4113231 0.1533231

A binomial model with known trials

The binomial branch keeps the same interface; the only addition is the known number of trials.

simb <- simulate_dynamic_binomial(n = 80, sigma = 0.12, trials = 50, seed = 4)
fit_bin <- fit_dynamic_model(simb$y, family = "binomial", trials = simb$trials,
                             horizon = 8, forecast_trials = 50,
                             nsave = NSAVE, nburn = NBURN, seed = 4)
summary(fit_bin)
#> Dynamic count model summary
#>   family = binomial | dynamics = rw | innovations = gaussian | zeros = none
#>   80 observations (0 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>            mean     sd   q2.5    q50  q97.5
#> innov_sd 0.2458 0.0317 0.1921 0.2442 0.3204
#> 
#> Fitted values: range of posterior means [24.73, 43.06]

plot_fitted(fit_bin)

Forecasting works the same way; supply forecast_trials (at fitting time) for the future trial sizes.

fc_bin <- forecast(fit_bin)
fc_bin$summary
#>   horizon   mean       sd q2.5 q50 q97.5
#> 1       1 42.009 3.558874   34  42    48
#> 2       2 41.941 3.814193   34  43    48
#> 3       3 42.092 3.873569   33  43    48
#> 4       4 42.041 4.128967   32  43    48
#> 5       5 42.190 4.302605   31  43    48
#> 6       6 42.263 4.589215   31  43    49
#> 7       7 42.224 4.764927   30  43    49
#> 8       8 41.942 5.152979   29  43    49

Zero inflation is available for the binomial family too: a structural-zero gate sits in front of the Binomial(m, p) process, exactly as for the Poisson. Set zero_inflation = TRUE (or zeros = "inflated") and read the per-zero diagnostics with structural_zero_prob().

simz <- simulate_dynamic_binomial(n = 80, sigma = 0.1, trials = 40, logit0 = 1.5,
                                  zero_inflation = 0.2, seed = 7)
fit_bz <- fit_dynamic_model(simz$y, family = "binomial", trials = 40,
                            zero_inflation = TRUE,
                            nsave = NSAVE, nburn = NBURN, seed = 7)
head(structural_zero_prob(fit_bz))
#>   time observed p_structural p_sampling
#> 1    3        0            1          0
#> 2    4        0            1          0
#> 3   14        0            1          0
#> 4   19        0            1          0
#> 5   24        0            1          0
#> 6   43        0            1          0

Choosing and changing priors

Every prior hyperparameter is exposed through dynamic_prior(); printing the object shows the current settings.

dynamic_prior()
#> <dynamic_prior>
#>   innovation variance ~ InvGamma(shape = 2.5, rate = 0.5)
#>   t degrees of freedom = 3 + Exp(mean = 6)
#>   mixture: 2 components, Dirichlet concentration = 1
#>   zero-inflation gate-open prob ~ Beta(1, 1)
#>   AR(1) rho ~ N(mean = 0, sd = 1) truncated to (-1, 1)  [ar1 only]
#>   drift/intercept mu ~ N(mean = 0, sd = 1)  [include_mu only]
#>   initial state ~ N(0, 100)  [rw and ar1]
#>   sv_prior: stochvol defaults

# a tighter prior favouring smoother latent paths
pr <- dynamic_prior(var_shape = 10, var_rate = 0.2)
fit_smooth <- fit_dynamic_model(sim$y, prior = pr,
                                nsave = NSAVE, nburn = NBURN, seed = 1)