Poisson Autoregression Models

• Poisson autoregression is a class of count time series models where the current count follows a Poisson distribution with intensities updated recursively from past data, covariates, or latent states.
• The framework includes observation-driven models such as PAR, INAR, and threshold autoregressions that capture persistence, regime changes, and non-linear dynamics in count data.
• Advanced variants extend to multivariate, copula-based, and time-varying formulations, enabling applications in fields like epidemiology, finance, and network analysis while addressing issues like overdispersion and limited dependency structures.

Searching arXiv for recent and foundational papers on Poisson autoregression and closely related count autoregressive models. arxiv_search(query="Poisson autoregression count time series INGARCH INAR threshold multivariate", max_results=10, sort_by="relevance") Searching for "Poisson autoregression" and related arXiv papers. Poisson autoregression denotes a family of count time-series models in which the conditional law of the observed process is Poisson and the conditional mean, or intensity, evolves recursively from past observations, past intensities, covariates, or latent states. In its standard observation-driven form, the model is written as

$$\mathcal{L}(Y_t \mid \mathcal{F}_{t-1})=\mathrm{Poisson}(\lambda_t), \qquad \lambda_t=\delta+\alpha \lambda_{t-1}+\beta Y_{t-1},$$

while broader usage includes Poisson autoregressive models with exogenous covariates, threshold and time-varying specifications, multivariate and network extensions, and thinning-based constructions such as INAR and INARMA that preserve Poisson marginals exactly (Wang et al., 2013). The literature also draws an important distinction between models that are conditionally Poisson and models whose unconditional or marginal distribution is exactly Poisson, because if $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$ and $\lambda_t$ is random, then the marginal distribution of $Y_t$ is typically not Poisson (Kong et al., 2023).

1. Canonical observation-driven formulations

The standard Poisson autoregressive framework is observation-driven: the current count is conditionally Poisson, and the conditional intensity is updated recursively from past data. In univariate PoARX form,

$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$

so past realized counts, past conditional means, and exogenous covariates all contribute to the intensity (Halliday et al., 2018). This is the count analogue of ARMA/GARCH-type dynamics, but with Poisson conditional distributions.

A related formulation appears in the PAR model of Brandt and Williams as used in the multiple-series estimation paper,

$$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$

which emphasizes mean reversion through the factor $1-\sum \rho_i$ and admits a convenient additive decomposition $m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$ (Redondo et al., 2021). In multivariate settings, the same basic logic is applied componentwise, either with separate intensity recursions joined by a copula or through fully vector-valued autoregressive dynamics (Doukhan et al., 2017).

This class is attractive because it handles discreteness through the Poisson measurement equation and serial dependence through the autoregressive mean. At the same time, the classical single-regime specification is limited. The SETPAR paper states explicitly that standard Poisson autoregression cannot model negative correlation in the observed counts, because the autoregressive effect is always of the same sign (Wang et al., 2013). The review of Poisson count time series adds a distinct limitation: conditional Poisson autoregression generally does not preserve Poisson marginals (Kong et al., 2023).

2. Poisson marginals, thinning, and integer-valued autoregression

A different branch of the literature seeks autoregressive dependence while retaining exact Poisson marginals. The fundamental construction is binomial thinning,

$$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$

which leads to the first-order integer autoregression

$X_t=\alpha\circ X_{t-1}+\epsilon_t,$

with $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$0 iid Poisson with mean $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$1. This yields a stationary Poisson marginal distribution with mean $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$2, and its autocorrelation is $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$3, so only nonnegative dependence is allowed (Kong et al., 2023).

The same review discusses CINAR models,

$Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$4

where a multinomial “decision vector” selects which lag contributes at time $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$5. These models preserve Poisson marginals but still do not generate negative autocorrelation (Kong et al., 2023). This makes thinning-based models the most literal count-autoregressive recursions with exact Poisson marginals, but also places clear restrictions on feasible dependence structures.

A more elaborate Poisson-marginal construction is the Poisson INARMA$Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$6 model,

$Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$7

with $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$8, $Y_t\mid \lambda_t\sim \text{Poisson}(\lambda_t)$9, and $\lambda_t$0. The model introduces an unobserved latent count process $\lambda_t$1, interpreted in the paper as an “INAR(1) with hidden juveniles,” and it has the equivalent distributional form

$\lambda_t$2

Its marginal mean and variance coincide,

$\lambda_t$3

and its autocorrelation function is

$\lambda_t$4

so it has ARMA$\lambda_t$5-type second-order structure rather than the simpler AR$\lambda_t$6-type correlation of INAR$\lambda_t$7 (Bracher, 2019).

The same paper shows that the process is equivalent to a binomially thinned INAR$\lambda_t$8 process through a latent Markov chain,

$\lambda_t$9

with $Y_t$0. This representation makes the mean, variance, correlation function, and time-reversibility straightforward, and it positions the model between INAR and INGARCH ideas: it preserves integer-valuedness and Poisson marginals while producing richer serial dependence (Bracher, 2019).

3. Thresholds, change points, and time variation

A major development beyond the classical single-recursion model is the introduction of regime-dependent intensity dynamics. The self-excited threshold Poisson autoregression (SETPAR) keeps the conditional Poisson observation law,

$Y_t$1

but replaces the single mean recursion with the two-regime update

$Y_t$2

with $Y_t$3, and $Y_t$4 (Wang et al., 2013). The paper’s motivation is explicit: the threshold structure remedies one of the drawbacks of classical Poisson autoregression by allowing possibly negative correlation in the observations. It also allows an “explosive” lower regime in the sense that $Y_t$5 may still be possible while the global process remains stable.

A different form of nonstationarity is piecewise parameter constancy. The piecewise autoregression paper models the conditional mean as

$Y_t$6

where $Y_t$7 is piecewise constant over time, with breakpoints $Y_t$8. Within each segment, the process follows an ordinary integer-valued autoregressive recursion; across segments, the parameter may jump (Diop et al., 2019). The framework explicitly covers Poisson-INARCH$Y_t$9, Poisson-INGARCH$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$0, negative-binomial INGARCH$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$1, binary INARCH$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$2, and INARCH$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$3. This suggests that “Poisson autoregression” is often used not as a single model, but as a base mechanism embedded in larger change-point structures.

Time variation can also be continuous rather than piecewise. The TV-PARX model specifies

$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$4

with fixed $$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$5 and $$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$6, but time-varying $$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$7 and $$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$8 updated by score-driven recursions (Angelini et al., 2022). The paper emphasizes that this flexibility is valuable when the data are affected by exogenous shocks, policy interventions, or crises, because persistence and covariate effects are not forced to remain constant over the entire sample.

4. Multivariate, copula-based, and network formulations

Multivariate Poisson autoregression can be constructed by combining univariate Poisson autoregressive margins with a dependence model across series. In the multivariate PoARX framework,

$$Y_t \mid \mathcal{F}_{t-1}\sim \text{Poisson}(\lambda_t), \qquad \lambda_t=\omega+\sum_{l=1}^p \alpha_l Y_{t-l}+\sum_{l=1}^q \beta_l \lambda_{t-l}+\eta\cdot x_{t-1},$$9

and the joint conditional distribution is built by a copula with Poisson marginals (Halliday et al., 2018). The paper uses Frank’s copula and stresses a structural feature of the model: dependence across series is introduced through the joint distribution, not through cross-lagged terms in the intensity recursion.

The multivariate count autoregression paper develops a broader copula-based construction using the continuous waiting-time representation of a Poisson process rather than imposing a copula directly on a vector of counts. Its linear specification is

$$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$0

and its log-linear specification is

$$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$1

componentwise (Doukhan et al., 2017). The log-linear form does not impose positivity constraints on the parameters and can accommodate negative feedback effects more flexibly.

Network structure provides a different route to multivariate dependence. In the Bayesian Poisson network autoregression mixture model (PNARM),

$$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$2

where $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$3 is a network autoregressive predictor, $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$4 is a node-specific baseline or offset term, and $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$5 is a latent cluster label (Hung et al., 2024). The model combines network-constrained dependence, latent clustering, and heterogeneous node dynamics. The paper’s Ireland COVID-19 application reports that a 5-component finite-mixture PNARM achieved mean absolute scaled error $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$6, compared with $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$7 for DDP PNARM, $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$8 for GAGNAR, $$Y_t \mid m_t \sim \text{Poisson}(m_t), \qquad m_t=\sum_{i=1}^p \rho_i Y_{t-i} +\left(1-\sum_{i=1}^p \rho_i\right)\exp(\delta_0 + X_t' \delta),$$9 for PNAR raw counts, and $1-\sum \rho_i$0 for PNAR population-adjusted counts.

The periodic multivariate Poisson autoregression extends the network setting further by allowing periodic, possibly infinite-memory kernels. For Type I periodicity,

$1-\sum \rho_i$1

and in the network version $1-\sum \rho_i$2, where $1-\sum \rho_i$3 is a normalized adjacency matrix (Khabou et al., 3 Apr 2025). The Rotavirus application in Berlin reports BIC values of $1-\sum \rho_i$4 for the proposed seasonal model and $1-\sum \rho_i$5 for PNAR$1-\sum \rho_i$6, and at horizon $1-\sum \rho_i$7 weeks the seasonal model significantly outperformed PNAR$1-\sum \rho_i$8 in 5 of 12 districts.

A complementary perspective is supplied by a rare-events scaling result: aggregates of interacting binary autoregressive processes can converge to a Poisson autoregression. In the interactive model studied in 2026,

$1-\sum \rho_i$9

with limiting recursion

$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$0

This provides a micro-foundation for the INGARCH$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$1-type model from a large system of interacting Bernoulli units under rare-events scaling (Bykhovskaya et al., 15 Apr 2026).

5. Estimation, asymptotics, and robust testing

Likelihood-based inference is central throughout the literature, but the computational route depends strongly on model structure. In SETPAR, the approximate log-likelihood is based on initialized recursions,

$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$2

and the paper proves strong consistency of the MLE and, when the threshold is treated as known, asymptotic normality

$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$3

under its stated assumptions (Wang et al., 2013). The threshold estimator is integer-valued, and the paper notes that $m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$4 eventually equals the true $m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$5.

For hidden-state count models, direct likelihood evaluation can be intractable. The INARMA$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$6 paper therefore treats the model as a hidden Markov model. Because the latent state space is countably infinite, the likelihood is computed after truncating the hidden state space at a large maximum value $m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$7, and the resulting log-likelihood is maximized numerically (Bracher, 2019). This representation is the computational route to parameter estimation.

In models with unknown breakpoints, estimation is often based on Poisson quasi-likelihood rather than a fully specified conditional law. The piecewise autoregression paper uses the segmentwise contrast

$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$8

defines a penalized criterion

$m_t=f_1(Y_{t-1},\dots,Y_{t-p})+f_2(X_t)$9

and implements the search by dynamic programming with $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$0 complexity (Diop et al., 2019). The penalty parameter is calibrated using the slope heuristic.

In multivariate PoARX, the preferred estimator is inference functions for margins (IFM). The margins are estimated separately from

$$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$1

and then the dependence parameter is estimated from the copula log-likelihood. The paper proves asymptotic normality for the marginal estimators, the dependence estimator, and the full parameter vector, with a Godambe or sandwich covariance form (Halliday et al., 2018).

Structural change testing raises a different inferential issue, especially under contamination. The robust testing paper studies Poisson autoregressive models with outliers and replaces the score-test machinery by a density power divergence criterion. The minimum density power divergence estimator is

$$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$2

and the resulting test statistic $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$3 converges under the null to

$$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$4

When $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$5, the method reduces to the score test of Kang and Song (2017); when $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$6, empirical sizes remain much closer to 5% in the simulated contaminated cases (Kang et al., 2019).

6. Scope, applications, and recurrent limitations

Applications span epidemiology, seismology, finance, trade counts, building traffic, and disease surveillance. SETPAR was applied to the number of major earthquakes in the world and performed better than standard PAR in AIC, in-sample MSE, and out-of-sample MSE, while BIC slightly favored the simpler PAR model (Wang et al., 2013). TV-PARX was applied to daily COVID-19 infections in Italy and to US corporate defaults, where the time-varying specification outperformed the constant-parameter alternative in log-likelihood and information criteria for the COVID-19 series and was particularly useful around crisis periods for defaults (Angelini et al., 2022). Multivariate PoARX was applied to counts of people entering and exiting a building, and the copula PoARX with covariates achieved the best out-of-sample test log score among the four compared models (Halliday et al., 2018). The additive multiple-series PAR paper used monthly counts of daily increases in four Asian stock indices and reported overall MAPE $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$7, rMSE $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$8, and MAD around $$\alpha \circ X=\sum_{i=1}^{X} B_i,\qquad B_i \overset{iid}{\sim}\text{Bernoulli}(\alpha),$$9 (Redondo et al., 2021).

Several limitations recur across the literature. First, conditional Poisson modeling does not imply Poisson marginals; the review paper makes this point directly through

$X_t=\alpha\circ X_{t-1}+\epsilon_t,$0

when $X_t=\alpha\circ X_{t-1}+\epsilon_t,$1 is random (Kong et al., 2023). Second, many exact Poisson-marginal constructions, notably INAR and CINAR, permit only nonnegative autocorrelation (Kong et al., 2023). Third, pure Poisson variance can be too restrictive in applied multivariate settings: randomized PIT histograms for PNARM showed U-shapes, indicating overdispersion relative to Poisson and suggesting that a negative binomial or other overdispersed count model may be more appropriate in future work (Hung et al., 2024). Fourth, near nonstationarity remains difficult: in the additive/backfitting PAR study, both the hybrid estimator and the extended Kalman filter underestimated parameters in nearly nonstationary models (Redondo et al., 2021).

The literature also expands Poisson autoregression beyond nonnegative counts and beyond directly observed autoregressive recursions. The generalized Poisson difference INGARCH models signed integer-valued data through

$X_t=\alpha\circ X_{t-1}+\epsilon_t,$2

with stationarity condition $X_t=\alpha\circ X_{t-1}+\epsilon_t,$3 (Carallo et al., 2020). The Poisson–Gamma Dynamical System instead places autoregression in a latent gamma state,

$X_t=\alpha\circ X_{t-1}+\epsilon_t,$4

thereby combining multivariate autoregression, gamma–Poisson overdispersion, and Bayesian nonparametric shrinkage (Schein et al., 2017).

Taken together, these developments show that Poisson autoregression is not a single recursion but a model class organized around Poisson conditional laws and autoregressive dependence. Within that class, the main axes of variation are whether Poissonity is imposed conditionally or marginally, whether dependence is single-regime or regime-dependent, whether dynamics are univariate or multivariate, and whether heterogeneity is handled by thresholds, time variation, latent states, copulas, or network structure.