Dynamic Models for Poisson, Binomial and Multinomial Time Series

Introduction

DynCount fits Bayesian state-space models to count 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 it is switched on with include_mu = TRUE. Under the random walk it acts as a drift, and under AR(1) it is an intercept that is always included. The observations are linked to the latent trajectory through one of three observation models:

• Poisson with a log link, so that \(y_t \sim \mathrm{Poisson}(e^{z_t})\);
• Binomial with a logit link, so that \(y_t \sim \mathrm{Binomial}(m_t,\, \mathrm{logit}^{-1}(z_t))\), where the trials \(m_t\) are known;
• Multinomial, for choice counts over \(K\) categories with known totals \(N_t\). Here each non-baseline category \(k\) has its own latent additive-log-ratio series \(z_{t,k} = \log(p_{t,k}/p_{t,b})\), which gives \(K - 1\) latent processes (see below).

An optional known offset \(o_t\) may be added to the linear predictor of all three observation models. It acts as a log-exposure for the Poisson mean, \(e^{o_t + z_t}\), as a shift of the binomial logit, and as a per-category shift of the multinomial log-ratios. It is a fixed, user-supplied input, not part of the latent process \(z_t\), and defaults to zero.

The distribution of the increments \(\varepsilon_t = z_t - \mu - \rho z_{t-1}\) is controlled by the innovations argument. It can be Gaussian ("gaussian", the default), Student-t ("t"), a finite scale mixture of normals ("mixture") or a stochastic volatility process ("sv", which requires the stochvol package). For the Poisson and binomial families, zeros can be handled by zero inflation with a time-constant gate-open probability (zeros = "inflated") or treated as missing values (zeros = "missing").

The model is estimated by Metropolis-within-Gibbs MCMC. The latent states are updated with adaptive random-walk Metropolis steps that use their Gaussian Markov random field full conditionals. The innovation parameters, \(\mu\) and \(\rho\) are drawn by Gibbs steps, with a Metropolis step for the Student-t degrees of freedom. Forecasts are obtained after fitting by forward simulation from the posterior draws.

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

Note that 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. The effective number of draws can be much smaller than nsave. For real analyses, use longer chains and the full series, and check convergence as shown in the section on convergence below.

Simulating data

The simulation helpers generate data with a known latent path, which is useful for checking recovery. simulate_dynamic_poisson() returns the counts y, the latent log-rate log_rate and the Poisson mean rate.

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

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.1528 0.0243 0.1082 0.1517 0.2049
#> 
#> Fitted values: range of posterior means [10.94, 93.31]

For this model the only global parameter is innov_sd, the standard deviation of the latent increments. Its true value in the simulation is 0.18. 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)

predict() summarises the in-sample fit. With the default type = "mean" it returns the posterior of the mean of \(y_t\), and with type = "response" it returns posterior predictive replicates of \(y_t\).

head(predict(fit)$summary)
#>   time observed     mean       sd     q2.5      q50    q97.5
#> 1    1       10 11.58507 2.344105 7.554892 11.42956 16.91972
#> 2    2       15 11.88528 2.122722 8.180216 11.73251 16.39812
#> 3    3       13 11.70898 2.002064 8.311704 11.46163 16.28356
#> 4    4        9 10.93698 1.614310 8.058833 10.84032 14.20784
#> 5    5        7 10.97079 1.680634 7.930762 10.94188 14.39572
#> 6    6       12 12.19498 1.775691 8.884296 12.24329 16.03465

The draws themselves are stored in fit$draws. Its main components are the latent states z (a draws x time matrix aligned with the observations), the increment variances sig2, the fitted means fitted, the replicates yrep and the parameter draws (innov_var, rho, mu and, depending on the model, nu, pi_open, mix_weight, sv_phi and others). The summary rows are derived from these draws. For example, innov_sd is the square root of innov_var, t_df summarises nu, ar1_rho summarises rho, drift_mu or intercept_mu summarise mu, and gate_open_prob summarises pi_open. The full layout is documented in ?fit_dynamic_model and ?summary.dynamic_fit.

Forecasting

Forecasts are obtained by forward simulation. For every stored posterior draw, forecast() propagates the latent path from the last in-sample state with the state equation, drawing the increments from the fitted innovation structure. It then draws a response from the observation model at each simulated state, so the intervals reflect parameter, state and innovation uncertainty. Future states carry no likelihood, so this gives exact draws from the posterior predictive distribution, and the horizon can be chosen after fitting.

fc <- forecast(fit, horizon = 8, seed = 1)
fc                                 # prints the forecast path
#> <dynamic_forecast> poisson, horizon 8
#>  horizon   mean       sd   q2.5 q50   q97.5
#>        1 60.900 14.68440 34.000  60  92.000
#>        2 62.308 18.27871 34.000  60 103.025
#>        3 63.368 21.93298 29.975  60 119.025
#>        4 63.852 24.56293 25.000  61 121.050
#>        5 64.182 25.88960 26.975  60 122.025
#>        6 64.804 28.14288 25.000  59 129.000
#>        7 65.320 29.99472 23.975  60 136.000
#>        8 66.302 31.62187 21.975  61 147.075
fc$final                           # the single 8-step-ahead forecast
#>   horizon   mean       sd   q2.5 q50   q97.5
#> 1       8 66.302 31.62187 21.975  61 147.075
plot_forecast(fit, horizon = 8, seed = 1)

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). Alternatively, fit_dynamic_model(..., horizon = H) simulates an H-step forecast right after sampling and stores it in the fit, and forecast(fit) without a horizon then returns the stored forecast. If a fit holds no stored forecast, forecast() needs a horizon and stops with an error otherwise.

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. AR(1) always includes an intercept. The package enables include_mu automatically, which gives the process a non-zero stationary mean \(\mu / (1 - \rho)\). The random walk corresponds to \(\rho = 1\), which lies outside the AR(1) parameter space, so the two specifications are separate models rather than nested ones.

# a genuinely stationary AR(1) log-rate with stationary mean 4, so mu = 4 * (1 - rho)
sim_ar <- simulate_dynamic_poisson(n = 150, sigma = 0.2, log_rate0 = 4,
                                   rho = 0.9, mu = 0.4, seed = 3)
# no need to set include_mu, because AR(1) enables the intercept automatically
fit_ar <- fit_dynamic_model(sim_ar$y, family = "poisson", latent_dynamics = "ar1",
                            nsave = NSAVE, nburn = NBURN, seed = 3)
summary(fit_ar)                    # reports the posteriors of ar1_rho and intercept_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.1668 0.0185 0.1321 0.1659 0.2055
#> ar1_rho      0.9102 0.0392 0.8278 0.9126 0.9824
#> intercept_mu 0.3529 0.1561 0.0673 0.3418 0.6827
#> 
#> Fitted values: range of posterior means [22.13, 119.22]

For a random walk with drift, keep the default dynamics and set include_mu = TRUE. The drift then appears as drift_mu in the summary.

Offset

For Poisson data with varying exposure, pass a known offset (a log-exposure term), and the mean becomes \(\exp(\text{offset}_t + z_t)\). When forecasting, supply the future exposures as forecast_offset, either one value per horizon or a single value that is recycled. If a model has an offset and no forecast_offset is given, forecast() warns and assumes an offset of zero.

expo  <- log(runif(120, 50, 200))  # known exposure, e.g. population at risk
sim_o <- simulate_dynamic_poisson(n = 120, sigma = 0.12, log_rate0 = -3.5,
                                  offset = expo, seed = 4)
fit_o <- fit_dynamic_model(sim_o$y, family = "poisson", offset = expo,
                           nsave = NSAVE, nburn = NBURN, seed = 4)
forecast(fit_o, horizon = 6, forecast_offset = log(120), seed = 4)$final
#>   horizon   mean       sd q2.5 q50  q97.5
#> 1       6 13.659 7.099516    4  12 31.025

Example data

The package ships two real weekly count series of irregular maritime crossings, which are loaded on first use. uk_weekly covers English Channel crossings from ISO week 2018-W01 to 2025-W11 (376 weeks), and med_weekly covers Mediterranean crossings from 2015-W40 to 2025-W11 (494 weeks). Both have the columns week (the ISO week label), count and date (the Monday of the week).

str(uk_weekly)
#> 'data.frame':    376 obs. of  3 variables:
#>  $ week : chr  "2018-W01" "2018-W02" "2018-W03" "2018-W04" ...
#>  $ count: int  0 0 0 0 7 0 0 0 0 0 ...
#>  $ date : Date, format: "2018-01-01" "2018-01-08" ...
plot(med_weekly$date, med_weekly$count, type = "h", xlab = "week",
     ylab = "crossings", main = "Weekly Mediterranean crossings")

Heavy-tailed and time-varying innovations

The Mediterranean series has large counts and few zeros. With Gaussian increments, the occasional large week-to-week jump would inflate the innovation variance for the whole series. The Student-t innovation makes the latent path robust to such jumps. For a fast build we use the most recent 120 weeks.

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.5913 0.1556 1.3331 1.5792  1.9482
#> t_df     6.3387 3.3456 3.1800 5.3718 16.5819
#> 
#> Fitted values: range of posterior means [1.87, 14109.57]

The posterior of the degrees-of-freedom parameter t_df indicates how heavy the increment tails are, with smaller values indicating heavier tails. Its prior is set with df_min and df_mean_excess in dynamic_prior().

A finite scale mixture of normals is a more flexible alternative. The number of components is set with dynamic_prior(mix_components = ...) and defaults to two. The components are exchangeable and are not identified individually, so summary() reports only innov_sd, the marginal standard deviation of the increments. The draws of the component weights and variances are stored in fit$draws$mix_weight and fit$draws$mix_var.

fit_mix <- fit_dynamic_model(med, family = "poisson", innovations = "mixture",
                             nsave = NSAVE, nburn = NBURN, seed = 2)
summary(fit_mix)
#> Dynamic count model summary
#>   family = poisson | dynamics = rw | innovations = mixture | zeros = none
#>   120 observations (1 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>            mean     sd  q2.5   q50  q97.5
#> innov_sd 1.5355 0.1277 1.303 1.527 1.8136
#> 
#> Fitted values: range of posterior means [2.76, 14126.63]

Stochastic volatility lets the increment variance change over time. The log-variance follows an AR(1) process whose level, persistence and volatility are reported as sv_mu, sv_phi and sv_sigma, and the per-increment variances are stored in fit$draws$sig2. This option requires the stochvol package.

fit_sv <- fit_dynamic_model(med, family = "poisson", innovations = "sv",
                            nsave = NSAVE, nburn = NBURN, seed = 2)
summary(fit_sv)
#> Dynamic count model summary
#>   family = poisson | dynamics = rw | innovations = sv | zeros = none
#>   120 observations (1 zeros), 1000 posterior draws
#> 
#> Global parameters (posterior summaries):
#>            mean     sd    q2.5    q50  q97.5
#> innov_sd 1.7272 0.2105  1.3921 1.7040 2.1899
#> sv_mu    0.2633 0.4821 -0.8404 0.2904 1.1359
#> sv_phi   0.7886 0.0928  0.5878 0.7974 0.9440
#> sv_sigma 0.8051 0.1565  0.5160 0.8056 1.1180
#> 
#> Fitted values: range of posterior means [0.92, 14116.60]
vol <- sqrt(apply(fit_sv$draws$sig2, 2, median))
plot(vol, type = "l", xlab = "week", ylab = "increment SD (posterior median)",
     main = "Time-varying innovation SD")

Checking convergence

MCMC output should be checked before it is interpreted. The effective sample size measures how many independent draws an autocorrelated chain is worth, and the coda package computes it directly from the stored draws.

ess <- function(x) round(unname(coda::effectiveSize(x)))
c(innov_sd = ess(sqrt(fit$draws$innov_var)),
  z_40     = ess(fit$draws$z[, 40]),
  t_df     = ess(fit_med$draws$nu))
#> innov_sd     z_40     t_df 
#>       63       96      122

With the short chains used here, several of these values are only a fraction of the 1000 kept draws. The innovation standard deviation of a smooth latent path is typically the slowest quantity to mix, so increase nsave (or thin) until the effective sample sizes of the quantities of interest are comfortably large. Trace plots, such as plot(sqrt(fit$draws$innov_var), type = "l"), and several chains with different seeds are useful further checks.

Zero inflation and structural zeros

uk_weekly (English Channel crossings) has many zeros in its early weeks. Turning on zero inflation lets the model separate structural zeros from sampling zeros. A structural zero arises when a latent gate switches the count off, whereas a sampling zero is produced by the Poisson process itself. The gate is drawn separately for every week, while the probability that it is open is a single parameter that is constant over time. We use the zero-heavy early window of the series here.

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.9036 0.0960 0.7366 0.8910 1.1111
#> gate_open_prob 0.6985 0.0493 0.6006 0.6988 0.7890
#> 
#> Fitted values: range of posterior means [0.09, 190.90]

In fit_dynamic_model(), zero_inflation = TRUE is shorthand for zeros = "inflated". The summary row gate_open_prob is the posterior of the gate-open probability \(\pi_{\text{open}}\), so one minus it is the probability of a structural zero. A simpler alternative is zeros = "missing", which treats all observed zeros as missing values.

structural_zero_prob() reports, for each observed zero, the posterior probability that it is structural. By default it returns only the zero observations, and zeros_only = FALSE returns one row per observation. plot_zero_inflation() shows these probabilities as a bar chart with one bar per observed zero.

sz <- structural_zero_prob(fit_zip)
head(sz, 10)
#>    time observed p_structural p_sampling
#> 1     1        0        0.484      0.516
#> 2     2        0        0.543      0.457
#> 3     3        0        0.563      0.437
#> 4     4        0        0.708      0.292
#> 5     6        0        0.671      0.329
#> 6     7        0        0.544      0.456
#> 7     8        0        0.439      0.561
#> 8     9        0        0.399      0.601
#> 9    10        0        0.381      0.619
#> 10   11        0        0.355      0.645
plot_zero_inflation(fit_zip)

In the resulting table, a 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). By contrast, a p_structural near 0 marks a zero that is consistent with a genuinely low rate.

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

A binomial model with known trials

The binomial branch keeps the same interface, and the only addition is the known number of trials. A single value is recycled over time.

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,
                             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.1425 0.0295 0.0927 0.1411 0.1999
#> 
#> Fitted values: range of posterior means [25.41, 42.83]
plot_fitted(fit_bin)

Forecasting works the same way. Supply the future trial sizes as forecast_trials, either one value per horizon or a single value that is recycled. If they are omitted, the last observed number of trials is used.

fc_bin <- forecast(fit_bin, horizon = 8, forecast_trials = 50, seed = 4)
fc_bin$summary
#>   horizon   mean       sd q2.5 q50 q97.5
#> 1       1 41.784 3.145033   35  42    47
#> 2       2 41.815 3.081277   35  42    47
#> 3       3 41.909 3.317789   35  42    48
#> 4       4 41.663 3.519064   34  42    47
#> 5       5 41.721 3.613061   34  42    48
#> 6       6 41.701 3.786018   33  42    48
#> 7       7 41.544 3.980441   33  42    48
#> 8       8 41.447 4.167560   32  42    48

Zero inflation is available for the binomial family too. Exactly as for the Poisson, a structural-zero gate sits in front of the Binomial(m, p) process. To use it, set zero_inflation = TRUE (or zeros = "inflated") and read the per-zero diagnostics with structural_zero_prob(). Note that the simulators use their zero_inflation argument differently. In simulate_dynamic_binomial() it is the probability of a structural zero, here 0.2, which corresponds to a gate-open probability of 0.8.

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)
summary(fit_bz)$params
#>                     mean         sd       q2.5       q50     q97.5
#> innov_sd       0.1481387 0.04418200 0.09874596 0.1375478 0.2849292
#> gate_open_prob 0.8037141 0.04375668 0.71258755 0.8058957 0.8815981
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

Multinomial choice counts

When each period yields counts over \(K\) mutually exclusive categories, family = "multinomial" models the category shares dynamically. One category \(b\) is the baseline – by default the one with the largest total count – and each of the remaining \(K - 1\) categories has its own latent additive-log-ratio (ALR) series \[ z_{t,k} = \log \frac{p_{t,k}}{p_{t,b}}, \qquad p_{t,k} = \frac{e^{z_{t,k}}}{1 + \sum_{j \ne b} e^{z_{t,j}}}, \qquad y_t \sim \mathrm{Multinomial}(N_t, p_t), \] where the row totals \(N_t\) are treated as known. Every ALR series follows the selected latent_dynamics and innovations, but the series share no parameters. Instead, each has its own innovation variance, its own \(\rho\) and \(\mu\), and its own copy of the prior. The series therefore interact only through the multinomial likelihood, and the sampler updates each series in turn with the other categories held at their current values. Zero inflation is not available for this family, and rows with a total of zero are treated as missing.

The simulation below uses category C as the baseline. We pass the simulated baseline to the fit, so that the fitted log-ratios are on the same scale as the simulated ones. Without it the fit would use the default baseline, the category with the largest total count, which here is B.

sim_m <- simulate_dynamic_multinomial(n = 80, sigma = c(0.15, 0.08), trials = 250,
                                      alr0 = c(-1, 0.3), baseline = 3,
                                      categories = c("A", "B", "C"), seed = 5)
head(sim_m$y)
#>       A   B   C
#> [1,] 39 132  79
#> [2,] 25 114 111
#> [3,] 48 118  84
#> [4,] 38 115  97
#> [5,] 28 118 104
#> [6,] 44 116  90
colSums(sim_m$y)
#>     A     B     C 
#>  3296 10081  6623
fit_m <- fit_dynamic_model(sim_m$y, family = "multinomial",
                           baseline = sim_m$baseline,
                           nsave = NSAVE, nburn = NBURN, seed = 5)
fit_m
#> <dynamic_fit>
#>   family      : multinomial
#>   categories  : 3 (A, B, C); baseline = C
#>   dynamics    : rw  (rho = 1)
#>   innovations : gaussian
#>   zeros       : none
#>   observations: 80  (zero-total rows: 0)
#>   draws kept  : 1000
summary(fit_m)$params
#>                   mean         sd       q2.5        q50     q97.5
#> innov_sd[A] 0.18505776 0.03156893 0.13277892 0.18227783 0.2603390
#> innov_sd[B] 0.07261848 0.01437368 0.05100505 0.07025055 0.1061444

The true innovation standard deviations are 0.15 for A and 0.08 for B. Posterior draws of multinomial fits carry a trailing category dimension. For example, fit_m$draws$fitted_prob is a draws x time x K array of shares, and the latent parameters (innov_var, rho, mu, …) are draws x (K - 1) matrices named by category. predict() and forecast() return long-format summaries with a category column, and the plot functions draw one panel per category. For forecasts, forecast_trials gives the future totals. If it is omitted, the last non-zero row total is used.

head(predict(fit_m, type = "prob")$summary)
#>   time category observed      mean         sd      q2.5       q50     q97.5
#> 1    1        A    0.156 0.1436402 0.01792291 0.1082552 0.1426827 0.1778244
#> 2    2        A    0.100 0.1333023 0.01411823 0.1085621 0.1320271 0.1618523
#> 3    3        A    0.192 0.1592094 0.01606813 0.1297601 0.1588734 0.1903972
#> 4    4        A    0.152 0.1498739 0.01551025 0.1209839 0.1487282 0.1815924
#> 5    5        A    0.112 0.1396657 0.01538785 0.1092454 0.1401441 0.1690480
#> 6    6        A    0.176 0.1578097 0.01537631 0.1288524 0.1583275 0.1902902
forecast(fit_m, horizon = 6, forecast_trials = 250, seed = 5)$final
#>   horizon category    mean       sd   q2.5 q50  q97.5
#> 1       6        A  32.117 15.42177  9.000  29  68.05
#> 2       6        B 134.882 17.32825 99.975 136 167.00
#> 3       6        C  83.001 12.78423 60.000  83 110.00
plot_fitted(fit_m)

plot_latent(fit_m, category = "A")

Two practical notes. First, the model is not invariant to the choice of baseline, because the dynamics are placed on the log-ratios relative to the baseline. If the baseline’s own share moves a lot, every ALR series inherits that movement. It is therefore best to choose a large category with a stable share (the baseline argument accepts a column name or index). Second, with \(K = 2\) and the second column as baseline the model is exactly the binomial model of the previous section.

Choosing and changing priors

Every prior hyperparameter is exposed through dynamic_prior(), and printing the object shows the current settings. The default prior on the innovation variance is \(\mathrm{InvGamma}(0.01, 0.01)\). It is weakly informative for increment standard deviations of about 0.1 and above, but it is not scale-free. For very smooth series, with increment standard deviations of a few hundredths, the results can be sensitive to this prior, and a sensitivity check with a smaller var_rate is advisable. A larger var_rate favours rougher latent paths.

dynamic_prior()
#> <dynamic_prior>
#>   innovation variance ~ InvGamma(shape = 0.01, rate = 0.01)
#>   t degrees of freedom = 3 + Exp(mean = 6)
#>   mixture: 2 components, Dirichlet concentration = 1,
#>            component variances ~ InvGamma(shape = 2.5, rate = 0.5)
#>   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
# an informative prior favouring rougher latent paths
pr <- dynamic_prior(var_shape = 2.5, var_rate = 0.5)
fit_inf <- fit_dynamic_model(sim$y, family = "poisson", prior = pr,
                             nsave = NSAVE, nburn = NBURN, seed = 1)
rbind(default = summary(fit)$params["innov_sd", ],
      informative = summary(fit_inf)$params["innov_sd", ])
#>                  mean         sd      q2.5       q50     q97.5
#> default     0.1527807 0.02432522 0.1081547 0.1516789 0.2048539
#> informative 0.2206108 0.02523688 0.1739750 0.2177266 0.2735754

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