R/DynCount-package.R

In DynCount: Bayesian Dynamic Models for Count Time Series

#' DynCount: Bayesian Dynamic Models for Count Time Series
#'
#' The \pkg{DynCount} package fits state-space models for non-Gaussian
#' time series. A latent trajectory \eqn{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 \eqn{z_t}.
#'
#' 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.
#'
#' @section Main entry points:
#' \describe{
#'   \item{[fit_dynamic_model()]}{Fit a Poisson, binomial or multinomial
#'         dynamic model (random walk or stationary AR(1)).}
#'   \item{[forecast()] / [predict()]}{Posterior predictive forecasts for any
#'         horizon (by forward simulation from the posterior draws), or the
#'         in-sample fitted values and replicates.}
#'   \item{[summary()]}{Posterior summaries of a fitted model.}
#'   \item{[simulate_dynamic_poisson()], [simulate_dynamic_binomial()],
#'         [simulate_dynamic_multinomial()]}{Simulate data.}
#'   \item{Plotting}{[plot_latent()], [plot_fitted()], [plot_forecast()],
#'         [plot_zero_inflation()].}
#' }
#'
#' @section Latent dynamics:
#' The `latent_dynamics` argument of [fit_dynamic_model()] selects the GMRF
#' state evolution \eqn{z_t = \mu + \rho z_{t-1} + \varepsilon_t}:
#' \describe{
#'   \item{`"rw"`}{A first-order random walk, i.e. \eqn{\rho = 1} fixed (the
#'         default).}
#'   \item{`"ar1"`}{A stationary AR(1) process, always with an intercept
#'         (`include_mu = TRUE`).}
#' }
#' Both share the precision \eqn{P = D_\rho^\top \mathrm{diag}(1/\sigma^2)
#' D_\rho}, where \eqn{D_\rho} is the generalised first-difference operator; it
#' is symmetric tridiagonal, so its bands are formed and the latent full
#' conditionals evaluated in \eqn{O(N)}. A scalar \eqn{\mu} (set
#' `include_mu = TRUE`) adds a \emph{drift} (RW) or \emph{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
#' \eqn{N(\code{init\_mean}, \code{init\_var})} prior on the first latent
#' state (default \eqn{N(0, 100)}), under both dynamics. 
#'
#' @section Multinomial (choice-count) series:
#' With `family = "multinomial"` the data are an \eqn{n \times K} matrix of
#' category counts with known row totals \eqn{N_t}. One category \eqn{b}
#' (by default the one with the largest total count) is the baseline, and
#' each of the other \eqn{K - 1} categories has its own latent
#' additive-log-ratio series \eqn{z_{t,k} = \log(p_{t,k} / p_{t,b})}, so
#' that \eqn{p_t} is the softmax of \eqn{(z_{t,\cdot}, 0)} and
#' \eqn{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 \emph{no} parameters. Each series has its own innovation
#' variance (and auxiliary innovation parameters), its own \eqn{\rho} and
#' \eqn{\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 \eqn{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 \eqn{N_t = 0} are uninformative and
#' handled as missing. Posterior draws carry a trailing category dimension,
#' see [fit_dynamic_model()].
#'
#' @section The four innovation structures:
#' The `innovations` argument controls the distribution of the latent
#' increments \eqn{\varepsilon_t = z_t - \mu - \rho z_{t-1}} (with
#' \eqn{\rho = 1} for the random walk and \eqn{\mu = 0} unless a drift/intercept
#' is included):
#' \describe{
#'   \item{`"gaussian"`}{\eqn{\varepsilon_t \sim N(0, \sigma^2)} with a single,
#'         constant variance.}
#'   \item{`"t"`}{A Student-t scale mixture: \eqn{\varepsilon_t \sim t_\nu(0,
#'         \sigma^2)}, robust to occasional large jumps. The degrees of
#'         freedom \eqn{\nu} are estimated from the data.}
#'   \item{`"mixture"`}{A finite scale mixture of normals with
#'         `mix_components` components (see [dynamic_prior()]), resulting in a
#'         flexible increment distribution.}
#'   \item{`"sv"`}{Stochastic volatility: \eqn{\log\sigma^2_t} follows an AR(1)
#'         process (delegated to the \pkg{stochvol} package).}
#' }
#'
#' @section Zero inflation:
#' Set `zeros = "inflated"` (or `zero_inflation = TRUE`) to fit zero
#' inflation with the Poisson or the binomial family: the observed count is
#' \eqn{y_t = v_t \tilde y_t}, where the gate
#' \eqn{v_t \sim \mathrm{Bernoulli}(\pi_{\mathrm{open}})} decides whether an
#' observed zero is \emph{structural} (gate closed) or a \emph{sampling} zero
#' produced by the Poisson/binomial process (gate open). The gate-open
#' probability \eqn{\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:
#' \emph{unconditional} draws (`fitted`, `yrep`), which include the gate and
#' are the quantities to compare with observed data (posterior predictive
#' checks), and \emph{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()].
#'
#' @section 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.
#'
#' @references
#' Zens, G. and Bijak, J. (2026). Dynamic Count Models with Flexible Innovation
#' Processes for Irregular Maritime Migration. \emph{The Annals of Applied
#' Statistics}, 20(2), 1671--1690. \doi{10.1214/26-AOAS2171}.
#'
#' @aliases DynCount DynCount-package
#' @keywords internal
"_PACKAGE"

## usethis namespace: start
#' @importFrom stats dbinom dnorm dpois plogis qlogis quantile rbeta
#'   rgamma rnorm rpois runif sd rbinom dgamma dexp rmultinom
#' @importFrom graphics abline lines points polygon legend barplot
#' @importFrom grDevices adjustcolor
#' @importFrom utils head
## usethis namespace: end
NULL