---
title: Poisson Autoregression Models
url: https://www.emergentmind.com/topics/poisson-autoregression
type: topic
---

# Poisson Autoregression Models

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 [1307.4626]. 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 [2310.10798].

## 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 [1806.04892]. 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)\) [2104.13520]. 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 [1704.02097].

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 [1307.4626]. The review of Poisson count time series adds a distinct limitation: conditional Poisson autoregression generally does not preserve Poisson marginals [2310.10798].

## 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 \(\epsilon_t\) iid Poisson with mean \(\lambda(1-\alpha)\). This yields a stationary Poisson marginal distribution with mean \(\lambda\), and its autocorrelation is \(\alpha^h\), so only nonnegative dependence is allowed [2310.10798].

The same review discusses CINAR models,
\[
X_t = D_{t,1}(\alpha\circ X_{t-1})+\cdots+D_{t,r}(\alpha\circ X_{t-r})+\epsilon_t,
\]
where a multinomial “decision vector” selects which lag contributes at time \(t\). These models preserve Poisson marginals but still do not generate negative autocorrelation [2310.10798]. 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\((1,1)\) model,
\[
X_t = \phi \circ S_t + I_t, \qquad
S_t = S_{t-1} - (X_{t-1} - I_{t-1}) + \kappa \circ X_{t-1},
\]
with \(I_t \stackrel{\text{iid}}{\sim} \text{Pois}(\tau)\), \(0<\phi\le 1\), and \(0<\kappa<1\). The model introduces an unobserved latent count process \(S_t\), interpreted in the paper as an “INAR(1) with hidden juveniles,” and it has the equivalent distributional form
\[
S_t \overset{d}{=} (1-\phi)\circ S_{t-1}+\kappa\circ X_{t-1}.
\]
Its marginal mean and variance coincide,
\[
\mathbb{E}(X_t)=\operatorname{Var}(X_t)=\frac{\tau}{1-\kappa},
\]
and its autocorrelation function is
\[
\rho(h)=\phi\kappa\xi^{h-1},\qquad \xi = 1-\phi(1-\kappa),
\]
so it has ARMA\((1,1)\)-type second-order structure rather than the simpler AR\((1)\)-type correlation of INAR\((1)\) [1910.07244].

The same paper shows that the process is equivalent to a binomially thinned INAR\((1)\) process through a latent Markov chain,
\[
Y_t = J_t + \xi\circ Y_{t-1}, \qquad
\tilde{Y}_t\mid Y_t \sim \text{Bin}\!\left(Y_t,\frac{\phi\kappa}{\xi}\right),
\]
with \(J_t\stackrel{\text{iid}}{\sim}\text{Pois}\!\left(\frac{\tau\xi}{\kappa}\right)\). 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 [1910.07244].

## 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,
\[
\mathcal{L}(Y_t \mid \mathcal{F}_{t-1}) = \textrm{Poisson}(\lambda_t),
\]
but replaces the single mean recursion with the two-regime update
\[
\lambda_t= \left\{ \begin{array}{cc}
d_1+a_{1}\lambda_{t-1}+b_{1}Y_{t-1},&Y_{t-1}\leq  r,\\
d_2+a_{2}\lambda_{t-1}+b_{2}Y_{t-1},&Y_{t-1} > r,
\end{array} \right.
\]
with \(d_i>0,\ a_i>0,\ b_i>0,\ i=1,2\), and \(r\in \mathbb{N}\) [1307.4626]. 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 \(a_1+b_1>1\) 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
\[
\lambda_t = (Y_t\mid \mathcal{F}_{t-1}) = f(Y_{t-1},Y_{t-2},\ldots;\theta_t^*),
\]
where \(\theta_t^*\) is piecewise constant over time, with breakpoints \(0<t_1^*<\cdots<t_{K^*-1}^*<n\). Within each segment, the process follows an ordinary integer-valued autoregressive recursion; across segments, the parameter may jump [1911.00989]. The framework explicitly covers Poisson-INARCH\((1)\), Poisson-INGARCH\((1,1)\), negative-binomial INGARCH\((1,1)\), binary INARCH\((1)\), and INARCH\((\infty)\). 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
\[
\log \lambda_{t+1} = \omega+\beta \log \lambda_t +\alpha_{t+1}(y_t-\lambda_t)\lambda_t^{-1} +\gamma_{t+1}'x_t +\psi'd_t,
\]
with fixed \((\omega,\beta)\) and \(\psi\), but time-varying \(\alpha_t\) and \(\gamma_t\) updated by score-driven recursions [2207.11003]. 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^j \mid \mathcal{F}_{t-1} \sim \text{Poisson}(\lambda_t^j), \qquad
\lambda_t^j = \omega^j + \sum_{l=1}^p \alpha_l^j Y_{t-l}^j + \sum_{l=1}^q \beta_l^j \lambda_{t-l}^j + \eta^j \cdot x_{t-1}^j,
\]
and the joint conditional distribution is built by a copula with Poisson marginals [1806.04892]. 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_{i,t}\mid {\cal F}^{Y}_{t-1},\lambda \sim \text{Poisson}(\lambda_{i,t}),\qquad
{\lambda}_t = {\bf d}+{\bf A}{\lambda}_{t-1}+{\bf B}{\bf Y}_{t-1},
\]
and its log-linear specification is
\[
\nu_t = {\bf d}+{\bf A}\nu_{t-1}+{\bf B}\log({\bf Y}_{t-1}+{\bf 1}_p),\qquad \nu_t\equiv \log\lambda_t
\]
componentwise [1704.02097]. 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_{i,t}\mid \mathcal{F}_{t-1}, Z_i, \theta \sim \text{Poisson}(\lambda_{i,t}), \qquad
\lambda_{i,t} = \theta_{1,Z_i} v_i + \theta_{2,Z_i} X_{i,t-1} + \theta_{3,Z_i} Y_{i,t-1},
\]
where \(X_{i,t-1}\) is a network autoregressive predictor, \(v_i\) is a node-specific baseline or offset term, and \(Z_i\) is a latent cluster label [2411.14265]. 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 \(0.46\), compared with \(0.52\) for DDP PNARM, \(0.65\) for GAGNAR, \(0.70\) for PNAR raw counts, and \(0.68\) 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,
\[
Y_t = N_t(\lambda_t),\qquad
\lambda_t = \psi\!\left(\mu_t+\sum_{k=1}^{t-1}\phi^{(t)}_{k}Y_{t-k}\right),
\]
and in the network version \(\phi_k^{(t)}=\alpha_k^{(t)}I_d+\beta_k^{(t)}W\), where \(W\) is a normalized adjacency matrix [2504.02649]. The Rotavirus application in Berlin reports BIC values of \(-16861.78\) for the proposed seasonal model and \(-16159.77\) for PNAR\((1)\), and at horizon \(h=13\) weeks the seasonal model significantly outperformed PNAR\((1)\) 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,
\[
X_t(N)=\sum_{i=1}^N y_{i,t}\Rightarrow X_t,\qquad
X_t\mid \mathcal F^X_{t-1}\sim \text{Pois}(\lambda_t),
\]
with limiting recursion
\[
\lambda_t=\bar c+\bar\gamma\,X_{t-1}+\beta\,\lambda_{t-1}.
\]
This provides a micro-foundation for the INGARCH\((1,1)\)-type model from a large system of interacting Bernoulli units under rare-events scaling [2604.14394].

## 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,
\[
\tilde{\ell}(\theta)=\sum_{t=1}^n \tilde{\ell}_t(\theta), \qquad
\tilde{\ell}_t=-\tilde{\lambda}_{t}+Y_t \log(\tilde{\lambda}_{t}),
\]
and the paper proves strong consistency of the MLE and, when the threshold is treated as known, asymptotic normality
\[
\sqrt{n}(\hat{\theta}-\theta_0)\xrightarrow{d} N(0, G^{-1})
\]
under its stated assumptions [1307.4626]. The threshold estimator is integer-valued, and the paper notes that \(\hat r\) eventually equals the true \(r\).

For hidden-state count models, direct likelihood evaluation can be intractable. The INARMA\((1,1)\) 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 \(Y^{\max}\), and the resulting log-likelihood is maximized numerically [1910.07244]. 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
\[
\widehat L_n(T,\theta)=\sum_{t\in T}\big(Y_t\log\widehat\lambda_t(\theta)-\widehat\lambda_t(\theta)\big),
\]
defines a penalized criterion
\[
\widetilde J_n(K,\underline t,\underline\theta) =\widehat J_n(K,\underline t,\underline\theta)+\kappa_n K,
\]
and implements the search by dynamic programming with \(\mathcal O(n^2)\) complexity [1911.00989]. 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
\[
l_j(\theta^j)=\sum_{t=1}^n \log \Pr(y_t^j \mid \mathcal{F}_{t-1};\theta^j),
\]
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 [1806.04892].

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
\[
\hat{\theta}_{\alpha,n} = \arg\min_{\theta \in \Theta}\sum_{t=1}^n l_{\alpha,t}(\theta),
\]
and the resulting test statistic \(T_n^\alpha\) converges under the null to
\[
\sup_{0\le s\le 1} \|B_d(s)\|^2.
\]
When \(\alpha=0\), the method reduces to the score test of Kang and Song (2017); when \(\alpha>0\), empirical sizes remain much closer to 5% in the simulated contaminated cases [1908.11466].

## 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 [1307.4626]. 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 [2207.11003]. 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 [1806.04892]. The additive multiple-series PAR paper used monthly counts of daily increases in four Asian stock indices and reported overall MAPE \(24.97\%\), rMSE \(3.21\), and MAD around \(2.71\) [2104.13520].

Several limitations recur across the literature. First, conditional Poisson modeling does not imply Poisson marginals; the review paper makes this point directly through
\[
\operatorname{Var}(Y_t)=E[\Lambda_t]+\operatorname{Var}(\Lambda_t)>E[\Lambda_t]=E[Y_t]
\]
when \(\lambda_t\) is random [2310.10798]. Second, many exact Poisson-marginal constructions, notably INAR and CINAR, permit only nonnegative autocorrelation [2310.10798]. 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 [2411.14265]. 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 [2104.13520].

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
\[
Z_t \mid \mathcal F_{t-1} \sim GPD(\tilde\mu_t,\tilde\sigma_t^2,\lambda), \qquad
\tilde\mu_t(1-\lambda)=\mu_t = \alpha_0+\sum_{i=1}^p \alpha_i Z_{t-i} + \sum_{j=1}^q \beta_j \mu_{t-j},
\]
with stationarity condition \(\sum_{i=1}^p \alpha_i+\sum_{j=1}^q \beta_j < 1\) [2002.04470]. The Poisson–Gamma Dynamical System instead places autoregression in a latent gamma state,
\[
y_v^{(t)} \sim \mathrm{Pois}\!\left(\delta^{(t)} \sum_{k=1}^K \phi_{vk}\,\theta_k^{(t)}\right), \qquad
\theta_k^{(t)} \sim \mathrm{Gam}\!\left(\tau_0 \sum_{k_2=1}^K \pi_{k k_2}\,\theta_{k_2}^{(t-1)},\, \tau_0\right),
\]
thereby combining multivariate autoregression, gamma–Poisson overdispersion, and Bayesian nonparametric shrinkage [1701.05573].

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.

Source: https://www.emergentmind.com/topics/poisson-autoregression