--- title: "Dynamic Models for Poisson and Binomial Time Series" author: "Gregor Zens" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Dynamic Models for Poisson and Binomial Time Series} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2 ) set.seed(1) library(DynCount) # Short MCMC runs keep the vignette fast to build; use longer runs in practice. NSAVE <- 1000L NBURN <- 1000L ``` ## 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](https://doi.org/10.1214/26-AOAS2171). Note: the MCMC runs below use short chains (`nsave` = `r NSAVE`, `nburn` = `r NBURN`) 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. ```{r simulate} sim <- simulate_dynamic_poisson(n = 80, sigma = 0.18, log_rate0 = 2.5, seed = 1) str(sim, max.level = 1) 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. ```{r fit-poisson} fit <- fit_dynamic_model(sim$y, family = "poisson", nsave = NSAVE, nburn = NBURN, seed = 1) fit summary(fit) ``` `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. ```{r plot-fitted} plot_fitted(fit) ``` ```{r plot-latent} 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. ```{r forecast} 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 fc$final # the single 8-step-ahead forecast 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)\). ```{r ar1} # 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 ``` ## 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. ```{r offset} 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 ``` ## 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. ```{r med-t} 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) ``` 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. ```{r fit-zip} data(uk_weekly) uk <- uk_weekly$count[1:130] mean(uk == 0) # many zeros fit_zip <- fit_dynamic_model(uk, family = "poisson", zero_inflation = TRUE, nsave = NSAVE, nburn = NBURN, seed = 3) summary(fit_zip) ``` `structural_zero_prob()` reports, for each observed zero, the posterior probability that it is structural. ```{r structural} sz <- structural_zero_prob(fit_zip, zeros_only = TRUE) head(sz, 10) 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. ```{r ppc-zip} # 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 ``` ## A binomial model with known trials The binomial branch keeps the same interface; the only addition is the known number of `trials`. ```{r binomial} 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) plot_fitted(fit_bin) ``` Forecasting works the same way; supply `forecast_trials` (at fitting time) for the future trial sizes. ```{r binomial-forecast} fc_bin <- forecast(fit_bin) fc_bin$summary ``` 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()`. ```{r binomial-zip} 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)) ``` ## Choosing and changing priors Every prior hyperparameter is exposed through `dynamic_prior()`; printing the object shows the current settings. ```{r priors} dynamic_prior() # 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) ``` ## 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](https://doi.org/10.1214/26-AOAS2171)