Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
71,934 characters · 11 sections · 30 citation commands
Generalized Autoregressive Multivariate Models: From Binary to Poisson
\address[Anna Bykhovskaya]{Duke University} \email{[email removed]}
\address[Nour Meddahi]{Toulouse School of Economics} \email{[email removed]}
Binary time series taking values in $\{0,1\}$ are central to modeling dynamic economic and financial phenomena in which outcomes represent discrete events or decisions rather than continuous quantities. In economics, they arise naturally in dynamic binary choice models, where agents repeatedly decide whether to take actions such as purchasing, investing, or entering a market rust1987optimal, heckman1981heterogeneity. In finance, binary indicators capture rare but economically important events, including extreme return days or bear markets nyberg2013predicting, corporate defaults and credit downgrades duffie2009frailty, and the exercise of binary options. Macroeconomic forecasting provides another key application, where the variable of interest is binary, such as the onset of a recession kauppi2008predicting.
Beyond individual processes, multivariate binary time series offer a natural framework for modeling interdependent or networked decisions, where outcomes for one unit depend on others. Such models are useful for analyzing financial contagion, correlated defaults, and evolving networks of economic linkages. Allowing the network structure to vary over time further enriches the analysis by capturing changing patterns of influence, information diffusion, or systemic risk propagation (cf., dynamic network formation of graham2016homophily). As we show, such cross-sectional interactions, when aggregated across a large panel under a rare-events scaling, provide a natural micro-foundation for Poisson autoregression, a leading model for count time series.
We develop a general framework for multivariate binary time series in which each binary observation is treated as a realization of a Bernoulli variable whose success probability is time-varying and depends on both past outcomes and past success probabilities, in the spirit of the GARCH-type dynamics of bollerslev1986generalized. For our model we establish, under general and potentially nonlinear dynamics for the success probabilities, the existence of a stationary solution and, under a contraction condition on the probability update function, its uniqueness. A distinctive feature of our formulation, compared to, for example, kedem2005regression and de2011dynamic, is that the success probabilities may depend on latent past probabilities and not only on observed outcomes. This additional layer allows for richer dynamics and feedback effects that cannot be captured when probabilities evolve solely as functions of observed realizations.
Our model fits into several important classes. When specialized to a linear form, our model corresponds to the generalized autoregressive score (GAS) framework of creal2013generalized. However, most of the existing results for GAS models blasques2014stationarity, blasques2022maximum do not apply in our setting. The reason is that Bernoulli variables, expressible as indicators that a uniform random shock falls below a success probability, generate discontinuities that violate the smoothness and Lipschitz assumptions crucial for the GAS theory. Our model is also a member of the class of positive-valued time series with exponential-family conditional distributions studied, in the univariate case, by aknouche2021count and, in the multivariate case with one lag of dependent variables, by lee2023modeling. Our analysis departs from those contributions by focusing specifically on cross-sectional interactions and their effect on the aggregate number of successes, a question those papers do not address.
The stochastic properties of related binary models have been studied by moysiadis2014binary and fokianos2017binary, and some of our structural results, particularly the existence and uniqueness of a stationary solution, are analogous in spirit to results obtained there, reflecting the robustness of the underlying probabilistic argument. The present paper departs from that work in two important respects. First, we consider an $N$-dimensional panel with cross-sectional interactions. Second, and more substantially, we target the rare-events asymptotic regime, which requires success probabilities to be close to zero and motivates working directly with $[0,1]^{N}$ rather than the reparametrizing probabilities via an inverse distribution function. The cited papers apply a smooth monotone transformation, such as the logistic or normal CDF, that maps probabilities to an unbounded domain. Under these transformations, small probabilities are mapped to values near $-\infty$, making analysis in this region problematic. Moreover, when combined additively, these highly negative transformed values decrease further, behaving more like products of probabilities than sums, which complicates interpretation.
Our central result provides a micro-level foundation for the integer-valued GARCH (INGARCH) model rydberg2000modelling,streett2000some,ferland2006integer, also known as Poisson autoregression fokianos2009poisson, one of the most widely used models for count time series. Under a rare-events scaling in which individual success probabilities vanish as the cross-sectional dimension $N$ grows, we show that the aggregate number of successes converges, in finite-dimensional distributions, to a Poisson autoregression.
The practical motivation is clear: many economically important phenomena, such as defaults in large credit portfolios, the exercise of binary options across many positions, or episodes of heightened systemic risk, are individually rare, yet their cross-sectional aggregate incidence is of primary interest. This contrasts with the original motivation of rydberg2000modelling, which arose from discretizing time in a Poisson process.
Our setting yields a richer family of Poisson autoregressions than the parametric specifications considered in fokianos2009poisson. In particular, our recursion for the Poisson intensity admits genuinely nonlinear functions of both the aggregate count and the limiting intensity. We also show that the result extends to network interaction structures: when the network is doubly stochastic and sufficiently dense, heterogeneous interaction weights wash out asymptotically, and the same Poisson autoregressive limit obtains. This illustrates a broader phenomenon in which network topology becomes irrelevant at the aggregate level once connectivity is sufficiently rich.
Finally, we propose maximum likelihood estimation (MLE) and establish consistency and asymptotic normality under a boundedness condition on the success probabilities, which rules out the rare-events regime. In that regime the Poisson MLE applied to the aggregate count is the natural substitute, and under structural assumptions such as homogeneity, the aggregation theory supplies an explicit mapping from the Poisson parameters back to those of the underlying binary system.
We illustrate the framework using daily returns on S$\&$P 100 constituents over 2005--2025, defining tail events as idiosyncratic returns falling below the $5$th percentile. The binary and Poisson parameter estimates are close, consistent with the aggregation theory, and the filtered Poisson intensity tracks the sum of individual success probabilities closely in sample. In out-of-sample forecasting, the heterogeneous binary model outperforms the Poisson-calibrated specification, indicating that the granularity of the binary likelihood remains empirically meaningful.
The remainder of the paper is organized as follows. Section (ref) introduces the model, provides motivating examples, and establishes the existence and uniqueness of a stationary solution. Section (ref) analyzes the model's aggregate behavior, focusing on the total number of successes within a given period. Section (ref) derives the consistency and asymptotic normality of the maximum likelihood estimator. Section (ref) extends the model by allowing for additional covariates. Section (ref) applies the proposed framework to S$\&$P 100 data. Finally, Section (ref) concludes. All proofs appear in the appendices.
Consider an $N$-dimensional binary vector $\mathbf y_t=(y_{1,t},\ldots,y_{N,t})^{\mathsf T}$, where $\mathsf T$ denotes transposition. Let $\mathcal I_t=\sigma(\mathbf y_\tau,\tau \leq t)$ denote the information available at time $t$, that is, the $\sigma$-field generated by $\{\mathbf y_t,\mathbf y_{t-1},\mathbf y_{t-2},\ldots\}$. The components of $\mathbf y_t$ are assumed to be independent conditional on $\mathcal I_{t-1}$. Given $\mathcal I_{t-1}$, each component $y_{i,t}$, $i=1,\ldots,N$, is assumed to be a Bernoulli random variable with success probability $p_{i,t}$,
The vector of success probabilities $\mathbf p_t=(p_{1,t},\ldots,p_{N,t})$ may depend on its own past values $\mathbf p_{\tau},\,\tau=t-1,\ldots,t-s$, as well as on past realizations $\mathbf y_{\tau},\,\tau=t-1,\ldots,t-q$. Formally,
where $s$ and $q$ are finite integers controlling the memory of the process, and $\mathbf g(\cdot)$ is a continuous mapping from $[0,1]^{Ns}\times\{0,1\}^{Nq}$ to $[0,1]^N$. The general form of $\mathbf g(\cdot)$ allows for nonlinearities and interactions, whereby $y_{i,t}$ depends not only on its own past but also on the past values and success probabilities of all $N$ variables. The dependence of $\mathbf p_t$ on both $\mathbf y_{t-1}$ and $\mathbf p_{t-1}$ mirrors the GARCH structure of bollerslev1986generalized. We, therefore, refer to the dynamics in (ref)--(ref) as the Generalized Autoregressive Binary Process (GAB).
Bernoulli distribution belongs to the one-parameter exponential family and, thus, GAB model fits in the class of models considered in aknouche2021count,lee2023modeling.
Let us introduce some useful examples of the function $\mathbf g(\cdot)$. We begin with a simpler setting where $N=1$, so that $\mathbf p_t=p_{1,t}\equiv p_t$ and $\mathbf y_t=y_{1,t}\equiv y_t$.
All of the above examples can be straightforwardly generalized to $N>1$, where each time series has its own $p_{i,t}$, independent of the others with $j\neq i$. More interesting cases arise when the dynamics of $p_{i,t}$ depend on the outcomes or probabilities of other series. A few examples are provided below.
All of the above multivariate examples can be readily extended to incorporate multiple lags of outcomes and probabilities, analogous to the transition from GAB($1,1$) to GAB($s,q$). They can also be generalized to allow for nonlinearities that aggregate $p_{i,t}$ and/or $y_{i,t}$ across $i$, for instance, in the spirit of the peer-effect interaction terms of bykhovskaya2022time, where a nonlinear function of neighbors' past outcomes enters the dynamics.
Figures (ref) and (ref) visualize multivariate GAB dynamics for $y_{i,t}$ and $p_{i,t}$: dots at row $i$ indicate realizations of successes for series $i$, while the continuous trajectory above each row represents the corresponding success probability $p_{i,t}$. Figure (ref) corresponds to the interactive model (ref) with $N=5$ and strictly positive coefficients $\gamma_i$. The figure reveals pronounced co-movement in the success probabilities induced by the common component $\frac{1}{N}\sum_{j=1}^N y_{j,t-1}$, while the series-specific term $\alpha_i y_{i,t-1}$ introduces heterogeneity across series.
Figure (ref) illustrates, using two series, how the opposite effect can arise, with the success probabilities moving in opposite directions. In particular, when the success probability $p_{j,t}$ increases, the negative term $-0.3 p_{j,t}$ becomes more pronounced,\footnote{The coefficients are chosen to ensure that $p_{i,t} \in [0,1]$.} which in turn reduces the success probability $p_{i,t+1}$ for $i \neq j$.
We next analyze the stationarity properties of the GAB model (ref)–(ref) in the asymptotic regime $T \to \infty$, with the number of time series $N$ held fixed.
Notice that no additional conditions are imposed, implying that all GAB models admit a stationary solution regardless of parameter values. The key features ensuring this property are the bounded support of all variables and the continuity of the probability function $\mathbf g$. Consequently, all examples from the preceding subsection, including the nonlinear Example 3 and the interactive Examples 6,7, possess stationary solutions. In the case of the logit GAB model (Example 3), note in particular that it is not necessary to impose the classical condition $|\beta|<1$.
Since all variables are bounded by construction (outcomes in $\{0,1\}$ and probabilities in $[0,1]$), the process has finite moments of all orders. Hence, strict stationarity implies weak stationarity, leading to the following corollary.
The law of iterated expectations and (ref) imply that $\mathbb{E}y_{i,t}=\mathbb{E}p_{i,t}$ and $\mathbb{E}_{t-1}y_{i,t}=p_{i,t},$ where the subscript $t-1$ denotes conditioning on the information available at time $t-1$, $\mathcal I_{t-1}$. Assuming $\sum\limits_{j=1}^{s}\alpha_j+\sum\limits_{j=1}^{q}\beta_j<1$, the linear GAB($s,q$) yields
and, assuming $\gamma+\beta<1$, its $N$-dimensional variation (ref) yields $$ \mu:=\mathbb{E}y_{i,t}=\mathbb{E}p_{i,t}=\frac{\omega}{1-\gamma-\beta}, $$ which matches the unconditional variance formula for classical GARCH($s,q$) and GARCH($1,1$).
In the interactive model (ref), assuming $\alpha_i+\beta_i<1$ for all $i$ and $\frac1{N}\sum_{i=1}^N\frac{\gamma_i}{1-\alpha_i-\beta_i}<1$, $$\mu_i:=\mathbb{E}y_{i,t}=\mathbb{E}p_{i,t}=\omega_i+ \alpha_i \mu_{i}+ \gamma_i \frac{1}{N}\sum_{j=1}^N \mu_{j}+\beta_i \mu_{i} =\frac{\omega_i}{1-\alpha_i-\beta_i}+\frac{\gamma_i}{N(1-\alpha_i-\beta_i)}\sum_{j=1}^N \mu_{j},$$ so that
and
In the network GAB$(1,1)$ model (ref), similar argument yields, in matrix form with $A:=\mathrm{diag}(\alpha_i+\beta_i)$ and $\Gamma:=\mathrm{diag}(\gamma_i)$,
Since $W$ is row-stochastic, the row sums of $A+\Gamma W$ equal $\alpha_i+\beta_i+\gamma_i\leq 1-\omega_i$, so for $\omega_i>0$ the matrix $I-A-\Gamma W$ is strictly row-diagonally dominant and hence invertible, giving $\mu=(I-A-\Gamma W)^{-1}\omega$. For general $W$ this solution does not simplify further. The interactive model (ref) corresponds to $W_{ij}=1/N$, the complete graph, in which case the uniform structure allows $\sum_i\mu_i$ to be solved explicitly, yielding the closed-form expressions (ref)--(ref).
The linear GAB($s,q$) model with $\sum\limits_{j=1}^{s}\alpha_j+\sum\limits_{j=1}^{q}\beta_j=1$, and hence $\omega=0$, admits multiple stationary solutions. In particular, both $y_{it}=p_{it}\equiv0$ and $y_{it}=p_{it}\equiv1$ constitute trivial stationary paths. This is in contrast to the IGARCH case, where the unit-root condition destroys weak stationarity; here, bounded support guarantees that strict stationarity is preserved even on the boundary, though at the cost of uniqueness.
Assumption (ref) requires that the Lipschitz coefficients of $\mathbf g(\cdot)$, when arranged into the companion matrix $\Phi$, yield a stable system in the sense that the spectral radius $\rho(\Phi)<1$. This is the direct analogue of the stationarity condition for a vector autoregression of order $\max(s,q)$. A simple sufficient condition (and equivalent for $N=1$) is the row-sum bound $\sum_{j=1}^{N}\sum_{\tau=1}^{s}\beta_{ij}^{\tau} +\sum_{j=1}^{N}\sum_{\tau=1}^{q}\alpha_{ij}^{\tau}\leq K<1$ for all $i$, which implies $\rho(\Phi)\leq\|\Phi\|_\infty\leq K<1$. For the network GAB$(1,1)$ model (ref), the companion matrix reduces to the single $N\times N$ block $A+\Gamma W$, so Assumption (ref) specializes to $\rho(A+\Gamma W)<1$, which also implies invertibility of $I-A-\Gamma W$ and hence the existence of the closed-form mean (ref).
Alternatively, the contraction assumption can be stated in an aggregate form:
Neither Assumption (ref) nor Assumption (ref) implies the other. For example, when $N=2$ and $s=q=1$, let $g_1(\mathbf p,\mathbf y)=\tfrac{2}{3}p_1$ and $g_2(\mathbf p,\mathbf y)=0$, then Assumption (ref) holds, but Assumption (ref) fails since $2/3>1/2$. Conversely, let $g_i(\mathbf p,\mathbf y)=f_i(p_1)+f_i(p_2)+f_i(y_1)+f_i(y_2)$, where $f_1(x)=\min\{0.2,0.4x\}$, $f_2(x)=\max\{0,0.4(x-0.5)\}$. Then $g_i\in[0,0.8]$ and each $f_i$ is Lipschitz with the coefficient $0.4$. Moreover, the intervals on which $f_1$ and $f_2$ are nonconstant do not overlap, so bound (ref) holds with $K = 0.4$ (and no smaller constant is possible). However, under Assumption (ref), we must have $\alpha_{ij} = \beta_{ij} = 0.4$, implying $\Phi = 0.8 (1,1)^{\mathsf T}(1,1)$ and hence $\rho(\Phi) = 1.6$.
While we could not locate Theorem (ref) in the literature, related results have appeared in moysiadis2014binary,fokianos2017binary,aknouche2021count,lee2023modeling. The nontrivial aspect in the binary setting arises from the discontinuous nature of $y_{i,t}$. A jump in $y_{i,t}$ from $0$ to $1$ can lead to markedly different values of future probabilities $\mathbf p_{t+1}$. In particular, given the information set at time $t-1$, since $y_{i,t}=I(u_{i,t}\leq p_{i,t})$, there is a discontinuity at $u_{i,t}=p_{i,t}$, which induces discontinuous and non-Lipschitz behavior in $\mathbf p_{t+1}$. The key to overcoming this difficulty lies in a suitable coupling construction.
In this section, we study the behavior of the total number of successes at a given time, $X_t(N)$, given by
While the previous section analyzed the limit of the system as $T \to \infty$ with $N$ fixed, we now investigate the complementary asymptotic regime in which $N \to \infty$. To initialize the dynamics, we consider an initial condition $$\{\mathbf p_0,\mathbf p_{-1},\ldots,\mathbf p_{1-\max\{s,q\}}\}.$$ For $t=0,\ldots,1-\max\{s,q\}$, we supplement this initial condition by setting $y_{i,t}=I(u_{i,t}\leq p_{i,t})$, where the variables $u_{i,t}$ are Bernoulli and independent of one another (and independent of all $u_{i,t}$ appearing in Eq. (ref)). The initial probabilities $\mathbf p_0,\mathbf p_{-1},\ldots,\mathbf p_{1-\max\{s,q\}}$ are allowed to be random. In this case, all errors $u_{i,t}$ are assumed to be independent of the initial condition.
Our focus is on the case of rare events, where each $p_{i,t}$ tends to zero as $N\to\infty$. Studying such rare events is important because many applications involve outcomes that occur infrequently but whose aggregate frequency carries significant information. Examples include default counts in large credit portfolios, the number of extreme returns across assets, or spikes in systemic risk indicators. In these settings aggregation not only smooths out idiosyncratic noise but also highlights the limiting distributional behavior of the system, providing a tractable framework for inference and risk assessment.
Figure (ref) illustrates the behavior of the system when $N$ is large and successes are rare. Although the total number of series is $N = 50$, the aggregate number of successes remains small throughout the sample. By contrast, Figure (ref) considers a setting with only $N = 5$ series but substantially higher success probabilities. Despite these differences in cross-sectional size and individual success rates, the total number of successes is of comparable magnitude across the two data-generating processes.
As we show below, aggregation naturally leads to a Poisson autoregression/INGARCH, providing a micro-level theoretical foundation for the applicability of such models. In the linear case, the conditional mean of the aggregate count follows a GARCH-type recursion, recovering exactly the model of ferland2006integer,fokianos2009poisson. More generally, our nonlinear Theorem (ref) shows that richer micro-level interactions survive aggregation in a precise sense: the limiting intensity $\lambda_t$ inherits a genuinely nonlinear structure, yielding a broader class of Poisson autoregressions than the parametric nonlinear specifications considered in fokianos2009poisson.
Theorem (ref) shows that different specifications for $\mathbf p_t$ can generate identical aggregate dynamics. For example, suppose that in the first model $\gamma_i = \tfrac{1}{3}$ for all $i$, while in the second model $\gamma_i = \tfrac{2}{3}$ for half of the units and $\gamma_i = 0$ for the remainder, with all other coefficients identical across the two models. In both cases, the average coefficient is $\bar{\gamma} = \tfrac{1}{3}$, implying identical aggregate dynamics. However, at the panel level, the two models differ substantially: in the second model, half of the probabilities evolve deterministically given $\mathbf p_0$.
There are many cases when the initial condition (ref) is satisfied. For example, if $p_{i,0}=\frac{r_i}{N}$, where $r_i$ are i.i.d. with finite mean. Then ${\frac1{N}\max_{i=1,\ldots,N}r_i\xrightarrow[N\to\infty]{p}0}$ and ${\sum_{i=1}^N p_{i,0}=\frac1{N}\sum_{i=1}^N r_{i}\xrightarrow[N\to\infty]{p}\mathbb{E}r_1}$. Moreover, it is satisfied by the stationary solution, as shown below.
Theorem (ref) can be extended to allow for multiple lags and various nonlinearities.
When $s=k=1$, $f_{\alpha,i}(y)=a_i y$, $f_{\gamma,i}\equiv0$, $f_{\gamma}(y,p)=y$, Theorem (ref) reduces to Theorem (ref). Alternatively, one can set $f_{\gamma}\equiv0$ and $f_{\gamma,i}(y,p)=\gamma_i y$. The linear multi-lag extension of Theorem (ref) can be obtained by setting $k=1$, $f_{\alpha,i}(y_1,\ldots, y_s)=\sum_{\tau=1}^{s}a_{i,\tau} y_{\tau}$, $f_{\gamma,i}(y,p)\equiv0$, $f_{\gamma}(y_1,\ldots, y_s,p_1,\ldots, p_s)=\sum_{\tau=1}^{s}(\delta_{\tau}^y y_{\tau}+\delta_{\tau}^p p_{\tau})$, where $\delta_{\tau}^y,\delta_{\tau}^p\in\mathbb{R}_+$. Alternatively, one can set $f_{\gamma}\equiv0$ and $f_{\gamma,i}(y_1,\ldots, y_s,p_1,\ldots, p_s)=\gamma_i \sum_{\tau=1}^{s}(\delta_{\tau}^y y_{\tau}+\delta_{\tau}^p p_{\tau})$. Both representations recover the same limiting intensity $\lambda_t$, because for linear functions the prefactor $1/N$ can be taken out of the argument.
The term $\frac1{N^{\kappa}}f_{\alpha,i}(y_{i,t-1},\ldots,y_{i,t-s})$ is a nonlinear and multi-lag generalization of $\frac{a_i}{N^{\kappa}}y_{i,t-1}$. For example, it allows for across-time interactions $y_{i,t-1}y_{i,t-2}$. This term disappears in the limit due to the prefactor $N^{-\kappa}$. Without the prefactor the term can occasionally take non-negligible values (e.g., $y_{i,t-1}=1$), violating the requirement of the Poisson approximation that all probabilities are small. The term $f_{\gamma,i}(\cdot)$ is a nonlinear and multi-lag generalization of the interactions, which are governed by the average number of successes. The averages are $o_p(1)$, leading to the average gradient of $f_{\gamma,i}$ in the limit as the outcome of the first order Taylor approximation. Thus, these nonlinearities asymptotically become linear. Finally, the term $\frac1{N}\boldsymbol\gamma_i^{\mathsf T}f_{\gamma}(\cdot)$ is the one which remains nonlinear in the limit, due to its arguments converging to non-degenerate random variables. Unlike the $f_{\alpha,i}(\cdot)$ and $f_{\gamma,i}(\cdot)$ terms, where individual heterogeneity enters through $i$-indexed coefficients that average out in the limit, $f_{\gamma}(\cdot)$ is common across all $i$: every coordinate $i$ is exposed to the same aggregate quantity $f_{\gamma}\left(\sum\limits_{j=1}^{N}y_{j,t-1},\ldots,\sum\limits_{j=1}^{N}y_{j,t-s}, \sum\limits_{j=1}^{N}p_{j,t-1},\ldots,\sum\limits_{j=1}^{N}p_{j,t-s}\right)$, with individual heterogeneity entering only through $\boldsymbol\gamma_i$.
A distinctive feature of the aggregation in Theorems (ref) and (ref) is that each coordinate $i$ interacts with the equal-weight cross-sectional aggregates $\sum\limits_{j=1}^{N}y_{j,t}$ and $\sum\limits_{j=1}^{N}p_{j,t}$, leaving no room for heterogeneous interaction weights. In some special cases, such as the network GAB model, this restriction can be relaxed, as illustrated in Theorem (ref) below.
A row-normalized network matrix $W$ is doubly stochastic if, for example, for each $i$, $d_i^{in}=d_i^{out}=d$, where $d_i^{in}=\sum\limits_{j=1}^{N} \mathcal I(W_{ji}\neq0)$, $d_i^{out}=\sum\limits_{j=1}^{N} \mathcal I(W_{ij}\neq0)$ are in- and out-degrees of a vertex $i$. That means, each individual has $d$ incoming and $d$ outgoing links, i.e. the network is $d$-regular. The doubly stochastic structure, combined with homogeneity of $\gamma_i=\gamma$, ensures that the network weights wash out in the limit. Thus, the limiting intensity equation (ref) is identical to that of Theorem (ref) with $\gamma_i=\gamma$.
The condition on the out-degree growing with $N$ ensures that terms $\sum_j W_{ij}y_{j,t-1}$ are of order $o_p(1)$, mimicking the scaling of the equal-weight average $\frac1{N}\sum\limits_{j=1}^{N}y_{j,t}$ and thereby ensuring the Poisson approximation applies.
We begin by discussing estimation when success probabilities are bounded away from $0$ and $1$, the regime where traditional maximum likelihood methods apply, such as moysiadis2014binary,fokianos2017binary,aknouche2021count,lee2023modeling. We then briefly discuss what can be done in rare-event scenarios.
Suppose that the success probability function $\mathbf g(\cdot)$ depends on an unknown parameter $\theta_0 \in \mathbb{R}^m$, such as the coefficients in the models discussed in Section (ref). The parameter $\theta_0$ can be estimated by maximum likelihood. Since the innovations $u_{i,t}$ are independent and uniformly distributed on $[0,1]$, the likelihood and log-likelihood functions conditional on the initial values $\mathbf y_0,\ldots,\mathbf y_{-q+1},\mathbf p_0,\ldots,\mathbf p_{-s+1}$ take the form
In empirical applications the first $q$ realizations of the binary process $\mathbf y_t$ can be used as initial values. The corresponding probabilities $\mathbf p_t$, however, are unobserved. A practical approach is, therefore, to initialize them using the sample means of $\mathbf y_t$. Lipschitz-type contraction in Assumptions (ref), (ref) ensures that asymptotically this initialization bias disappears.\footnote{Let $\mathbf p_t$ and $\mathbf p'_t$ be two trajectories with different initial probability values (but with the same sequence of outcomes $\mathbf y$). Let $\ell_{i,t}(p):=y_{i,t}\log p+(1-y_{i,t})\log(1-p)$. Since the trajectory of $\mathbf y_t$ is the same, both Assumptions (ref) and (ref) imply that for all $i,t$, $|p_{i,t}-p'_{i,t}|\leq C\rho^t$ for some $C>0,\rho\in[0,1)$. Given that $g_i(\cdot;\theta)\in[\varepsilon,1-\varepsilon]$, as in Assumption (ref)(ref), $|\ell_{i,t}(p_{i,t})-\ell_{i,t}(p'_{i,t})|\leq\frac1{\varepsilon}|p_{i,t}-p'_{i,t}|\leq \frac{C}{\varepsilon}\rho^t$ and $\frac1{T}\sum_{i,t}|\ell_{i,t}(p_{i,t})-\ell_{i,t}(p'_{i,t})|\leq\frac{CN}{T\varepsilon}\frac{1}{1-\rho}\xrightarrow[T\to\infty]{}0$ uniformly in $\theta$ belonging to a compact.}
Continuity of $\mathbf g$ in $\theta$ and its mapping into $[0,1]^N$ imply that the parameter space $\Theta$ is closed. Boundedness, and hence compactness, is trivially satisfied in all linear GAB examples. For general nonlinear $\mathbf g$, compact support can be relaxed if, with positive probability, $\mathbf g\left(\{\mathbf p_{\tau}\}_{\tau=t-1,\ldots,t-s},\,\{\mathbf y_{\tau}\}_{\tau=t-1,\ldots,t-q};\theta\right)$ goes to $0$ or $1$ as $\theta \to \infty$. This ensures that the log-likelihood is never maximized at unbounded values of $\theta$.
Assumption (ref)(ref) is required to ensure that the log-likelihood is well-defined. It resembles the standard GARCH condition that the intercept be bounded away from zero (see, e.g., francq2019garch). In our setup, when the functions $g_i$ are linear, this means that all intercepts are separated from zero and that the sums of all coefficients in $g_i$ are separated from one. Note that if, for all $i$, $g_i(0,\ldots,0;\theta_0)=0$, then $y_{i,t}=p_{i,t}=0$ is a degenerate stationary solution. Similarly, if, for all $i$, $g_i(1,\ldots,1;\theta_0)=1$, then $y_{i,t}=p_{i,t}=1$ is a degenerate stationary solution. In both cases, the log-likelihood becomes infinite. By imposing (ref)(ref), we rule out such degenerate settings.
Assumption (ref)(ref) requires the existence of a unique parameter value $\theta$ that is compatible with the data, i.e., identification of $\theta_0$. This condition is again satisfied for linear models, since the stationary distribution of $\mathbf p_t$ is continuous (when it is not degenerate).
The maximum likelihood estimator solves $$\hat{\theta}^{MLE}=\arg\max\limits_{\theta\in\Theta} \mathcal{Q}_T(\theta).$$
Assumption (ref) or (ref) is imposed to ensure ergodicity and uniqueness of a stationary solution.
Note that the boundedness condition in Assumption (ref)(ref) rules out rare events, where $p_{i,t}$ is close to zero. Both consistency and asymptotic normality rely on this condition: it ensures the log-likelihood is well-defined and uniformly dominated, and keeps the score and Hessian bounded. As we illustrated in Section (ref), when $N$ is large and $p_{i,t}$ is small, the aggregate count $\sum_{i=1}^N y_{i,t}$ is approximately Poisson with mean $\sum_{i=1}^N p_{i,t}$. Thus, a natural alternative for estimation in such settings is the Poisson autoregressive model, which replaces the binary observation equation with a Poisson likelihood.
In the framework of Theorem (ref), when $\beta<1$, $\gamma_i=\nu c_i$ for all $i$ (that is, when some form of homogeneity is present), the following intuition can be used to identify the Poisson log-likelihood in the MLE maximization in (ref). First, in line with $\max_p p_{i,t}\xrightarrow[T\to\infty]{p}0$, we expect $$\left|\sum_{i}y_{i,t}\log(1-g_{i,t,\theta})\right|\leq|\log(1-\max_i g_{i,t,\theta})|\sum_{i}y_{i,t}\xrightarrow[T\to\infty]{p}0.$$ Therefore,
so that the first two terms operate on aggregate quantities and exactly match the Poisson log-likelihood. The last term is close to being independent of aggregate quantities and can therefore be maximized separately. To see this, applying iterative back-substitution yields $$p_{i,t}=\frac{c_i}{N(1-\beta)}+\beta^t\left(p_{i,0}-\frac{c_i}{N(1-\beta)}\right)+\frac{a_i}{N^{\kappa}}\sum\limits_{\tau=0}^{t-1}\beta^{\tau}y_{i,t-1-\tau} +\frac{\nu c_i}{N}\sum\limits_{\tau=0}^{t-1}\beta^{\tau}\sum\limits_{j=1}^Ny_{j,t-1-\tau},$$
When $t$ is large, $\beta^t\approx 0$, and the second terms in both expressions are negligible. Since $\kappa>0$, the third terms are also negligible, leaving $$\frac{g_{i,t,\theta}}{\sum\limits_{j=1}^{N}g_{j,t,\theta}} \approx\frac{\frac{c_i}{N}\left(\frac1{1-\beta}+\nu\sum\limits_{\tau=0}^{t-1}\beta^{\tau}\sum\limits_{j=1}^Ny_{j,t-1-\tau}\right)} {\frac{\sum_j c_j}{N}\left(\frac1{1-\beta}+\nu\sum\limits_{\tau=0}^{t-1}\beta^{\tau}\sum\limits_{j=1}^Ny_{j,t-1-\tau}\right)}=\frac{c_i}{\sum_j c_j}.$$ That is, only the relative parameters $\tilde{c}_i:=c_i/\sum_j c_j$ remain. These lie on a simplex and are independent of the aggregate values, so the original maximization problem splits into two separate problems, the first of which corresponds to a Poisson likelihood.
The model can be extended to include exogenous covariates $\mathbf z_t\in\mathbb{R}^k$, so that $$p_{i,t}=g_i\left(\{\mathbf p_{\tau}\}_{\tau=t-1,\ldots,t-s},\,\{\mathbf y_{\tau}\}_{\tau=t-1,\ldots,t-q},\,\mathbf z_t\right).$$ Provided $\mathbf z_t$ is strictly stationary and geometrically ergodic, Theorems (ref) and (ref) hold unchanged under Assumptions (ref)--(ref) as stated. This follows because in the coupling construction both trajectories share the same realization of $\mathbf z_t$, so it cancels identically from the Lipschitz bounds (ref) and (ref). Similarly, estimation Theorems (ref) and (ref) continue to hold.
The more interesting question is how to incorporate covariates in the aggregation framework. If $k$ is fixed, then $\mathbf z_t$ must enter each equation for $p_{i,t}$ with a coefficient vanishing in $N$. Otherwise, the finite order of $\mathbf z_t$ would prevent $p_{i,t}$ from being of order $1/N$. Paralleling the treatment of the intercept $\omega_i$, we therefore add a linear term $\frac{1}{N}\delta_i^{\mathsf T}\mathbf z_{t}$ to each $p_{i,t}$ equation in interactive and network GAB models. Letting $\bar{\delta}:=\lim_{N\to\infty}\frac1{N}\sum_{i=1}^{N}\delta_i$, this contributes an extra term $\bar{\delta}^{\mathsf T}\mathbf z_t$ to the right-hand side of the evolution equation for $\lambda_t$.
If $k$ grows with $N$, the covariates must themselves aggregate to a limiting quantity. We decompose $\mathbf z_t=(\mathbf z_t^{com},\mathbf z_t^{ind})$, where $\mathbf z_t^{com}\in\mathbb{R}^{r_c}$ is a common component of fixed dimension affecting all units, and $\mathbf z_t^{ind}=(\mathbf z_{1,t}^{ind},\ldots,\mathbf z_{N,t}^{ind})\in\mathbb{R}^{Nr}$ stacks unit-specific components with $\mathbf z_{i,t}^{ind}\in\mathbb{R}^r$, so that $k=Nr+r_c$. The common component is handled exactly as in the fixed $k$ case. For the individual components, $p_{i,t}$ may depend on $\mathbf z_{i,t}^{ind}$ directly and on the cross-sectional average $\frac{1}{N}\sum_{j=1}^N\mathbf z_{j,t}^{ind}$ (or the weighted average $\sum_{j=1}^N W_{ij}\mathbf z_{j,t}^{ind}$ in the network GAB case, as in Theorem (ref)). The direct effect of $\mathbf z_{i,t}^{ind}$ on $p_{i,t}$ must be scaled by $1/N$, for the same order-$1/N$ reason as above, and must appear with a homogeneous coefficient $\eta$, so that the aggregation yields $\eta\frac1{N}\sum_{i=1}^N\mathbf z_{i,t}^{ind}$. The cross-sectional average must carry an additional factor of $1/N$, since unlike $y_{i,t}$ and $p_{i,t}$, the scale of the exogenous $\mathbf z_{i,t}^{ind}$ is not reduced by aggregation; the corresponding term in $p_{i,t}$ is therefore $\frac{\zeta_i}{N}\frac{1}{N}\sum_{j=1}^N\mathbf z_{j,t}^{ind}$ (or nonlinear bounded\footnote{Boundedness guarantees that $p_{i,t}\in[0,1]$ and that $\mathbb{E}f_{z}\left(\frac{1}{N}\sum_{j=1}^N\mathbf z_{j,t}^{ind}\right)\xrightarrow[N\to\infty]{}\mathbb{E}f_{z}\left(\bar{\mathbf z_{t}}\right)$.} function $\frac{\zeta_i}{N}f_{z}\left(\frac{1}{N}\sum_{j=1}^N\mathbf z_{j,t}^{ind}\right)$, so that $f_z$ survives in the limit analogously to the $f_{\gamma}$ term in Theorem (ref)). Thus, assuming that $\frac{1}{N}\sum_{i=1}^N\mathbf z_{i,t}^{ind}\xrightarrow[N\to\infty]{p}\bar{\mathbf z}_t$, $\mathbb{E}\frac{1}{N}\sum_{i=1}^N\mathbf z_{i,t}^{ind}\xrightarrow[N\to\infty]{}\mathbb{E}\bar{\mathbf z}_t<\infty$, and $\max_i\|\mathbf z_{i,t}^{ind}\|/N\xrightarrow[N\to\infty]{p}0$ for each $t$, the aggregated covariate $\bar{\mathbf z}_t$ appears on the right-hand side of the limiting evolution equation for $\lambda_t$.
In this section, we apply the interactive GAB model to study periods of unusually low stock returns among the constituents of the S$\&$P 100. Such episodes tend to coincide with market stress and elevated financial risk. We use binary indicators of tail events to study how downside risks emerge, cluster over time, and propagate across assets, defining a tail event for stock $i$ as its return falling below the $5$th percentile of its empirical distribution.
Our sample consists of daily returns for $87$ S$\&$P 100 constituents available over the full period 01.01.2005--01.01.2025, yielding $5,032$ time series observations. The $87$ stocks are those for which data are available throughout the entire span; we also retain only one of the two Google share classes. Further details on the sample composition are provided in Section (ref) of the Appendix.
We subtract the risk-free rate and regress returns on the Fama-French five factors, then work with the idiosyncratic residuals. The final year of data is reserved for forecast evaluation, so estimation uses observations up to 01.01.2024, giving $T= 4,780$. The $5$th percentile threshold for each stock is computed on this estimation sample, and $y_{i,t} = 1$ whenever stock $i$'s idiosyncratic return falls below its threshold, and $y_{i,t} = 0$ otherwise. Figure (ref) summarizes the resulting binary panel. As expected, tail events cluster around the crisis periods: financial crisis, COVID, and Fed's monetary policy tightening. This clustering motivates the inclusion of the cross-sectional interaction term as in Eq. (ref), which captures the market-wide stress and allows individual tail probabilities to adjust.
We first estimate the interactive GAB model (ref). The resulting estimates are shown in Figure (ref). In line with the aggregation framework, both $\omega_i$ and $\alpha_i$ are small, with $\text{mean}(\omega_i)=0.0006,\,\text{mean}(\alpha_i)=0.027$. The averages of the remaining coefficients are $\text{mean}(\gamma_i)=0.18,\, \text{mean}(\beta_i)=0.78$.
We then re-estimate the model imposing $\alpha_i=0$, since this coefficient vanishes under aggregation, and compare the results with the Poisson MLE based on the aggregated data. The results are summarized in Table (ref). They are broadly consistent with the aggregation theorems, with the Poisson coefficients close to the averages of the GAB estimates. When $\alpha_i=0$ is imposed, the intercept absorbs the mean contribution previously attributed to $\alpha_i$, leading to $\sum_i\omega_i$ larger than $\bar{c}$. The overall similarity is further supported by the close agreement between the filtered Poisson intensity $\lambda_t$ and the sum of GAB probabilities $\sum_i p_{i,t}$, as shown in Figure (ref).
Finally, we compare the forecasting performance of the GAB model using the one-year sample that was excluded from the estimation stage. The out-of-sample binary data are constructed by comparing returns to the same $5$th percentile threshold used previously. We generate recursive 1-step-ahead forecasts of $y_{i,t}$ over the out-of-sample period ($T+1,\ldots,T+252$, one year of trading days) based on four models.
The first model is the interactive GAB specification without $\alpha_i$ (heterogeneous parameters). We impose $\alpha_i=0$ to align the GAB specification with the Poisson limit structure. The second is an interactive GAB model with homogeneous parameters calibrated from the Poisson maximum likelihood estimates, i.e., $\omega_i = \bar{c}/N$, $\gamma_i = \bar{\gamma}$, and $\beta_i = \beta$. The third is a constant forecast $\hat{y}_{i,t} = 0.05$, reflecting the unconditional $5\%$ probability of observing a success. The fourth is a persistence-based forecast $\hat{y}_{i,t} = y_{i,t-1}$.
We report the mean squared errors of these forecasts in Table (ref). The results show that the first model, the interactive GAB with heterogeneous parameters, achieves the best forecasting performance. Although the GAB and Poisson estimation results in Table (ref) are close, the second (homogeneous) model performs slightly worse. Both naive benchmarks are substantially outperformed.
Figure (ref) presents the out-of-sample data, both aggregated and disaggregated, together with the path of $\sum_i p_{i,t}$ for the first two models. For the second model, this path coincides with the trajectory of the Poisson intensity $\lambda_t$. The out-of-sample period features a large number of abnormally low returns ($21.3\%$ of the data points\footnote{This is consistent with MSE results. For the fourth model, the squared forecast error is equal to $0$ when two consecutive realizations coincide and $1$ when they differ. If the observations were i.i.d. Bernoulli, the latter would occur with probability $2p(1-p) \approx 0.32$ for $p \approx 0.2$. For the third model, the squared forecast error takes values $0.05^2$ and $0.95^2$, and $MSE\approx 0.2\cdot0.05^2+0.8\cdot0.95^2\approx0.18$. For the first two models, $MSE\approx0.16$ is close to the optimal value, which corresponds to the forecast $\hat{y}_{i,t}\approx0.2$ yielding $MSE=p(1-p)^2+(1-p)p^2\approx0.16$.}), indicating elevated uncertainty, likely reflecting elevated macroeconomic uncertainty in 2024. This highlights the importance of the interaction term $\frac{1}{N}\sum_i y_{i,t}$ in capturing overall market conditions; the term which is absent in alternative third and fourth models.
This paper develops a general class of multivariate binary autoregressive models, in which each observation is a Bernoulli variable whose success probability evolves according to a GARCH-type recursion driven by past outcomes and past probabilities. The framework accommodates nonlinear dynamics, cross-sectional interactions, and network-weighted dependence structures.
Three sets of results are established. The bounded support of binary outcomes and probabilities guarantees, without any parameter restrictions, the existence of a strictly stationary solution, a feature with no counterpart in continuous-valued GARCH models, while uniqueness and geometric ergodicity follow from a coupling argument that exploits a Lipschitz contraction condition on the probability update function. The most distinctive results concern aggregation: under a rare-events scaling in which individual success probabilities are of order $1/N$, the aggregate number of successes converges to a Poisson autoregression, providing a rigorous micro-foundation for this commonly used model. The network extension shows further that, under sufficient connectivity, heterogeneous interaction weights wash out asymptotically, leaving a limiting intensity equation that does not depend on the network topology. Finally, in the moderate-probability regime, the MLE is consistent and asymptotically normal with asymptotic variance equal to the inverse Fisher information; in the rare-events regime, the Poisson likelihood applied to the aggregate count is the natural substitute, and the aggregation theory supplies the mapping from its parameters back to those of the underlying binary system.
Several directions remain open. On the theoretical side, a unified asymptotic theory in which both $T$ and $N$ grow simultaneously would constitute a valuable generalization. The extension of asymptotic normality for the MLE to settings where $N$ grows with $T$ is a natural complement. At the intersection of the two regimes, a composite likelihood that exploits both the individual binary structure and the aggregate Poisson approximation may yield efficiency gains relative to either likelihood alone. On the empirical side, the binary panel framework is well suited to modeling correlated binary decisions in large administrative datasets: voting records, loan defaults, and insurance claims. The network GAB model in particular provides a tractable vehicle for studying the role of connectivity in propagating financial or economic shocks.