Poisson counts with a switching rate

Events whose arrival rate follows a hidden regime, and the distribution of their count.

The model

Events arrive one at a time, as insurance claims, defaults in a portfolio or requests to a server do.

Their rate depends on a hidden regime $y_t \in \{1, 2\}$: while $y_t = i$ they arrive as a Poisson process with rate $\ell_i$.

The regime switches at rate $\lambda$ in each direction. This is a Markov-modulated Poisson process.

There are two time scales: events arrive every $1/\ell_i$ on average, and the regime changes every $1/\lambda$.

Parameters: $\ell = (8, 1)$ per unit time, the process starts in the busy regime, and $N_t$ counts the events up to $t = 1$.

The quantity

The whole distribution $\Pr(N_t = k)$, $k = 0, 1, 2, \dots$, through the generating function

Averaging

When the regime switches infinitely fast, events arrive as a Poisson process with the average rate $\bar\ell = 4.5$, and the variance of the count equals its mean. Switching at a finite rate makes the count more variable than Poisson, and the corrector carries that overdispersion.

Reduction to a linear system

Condition on what happens in a short interval $dt$ at the start. While the regime is $i$, an event occurs with probability $\ell_i\,dt$ and multiplies $z^{N}$ by $z$; the regime jumps to $j$ with probability $Q_{ij}\,dt$. So

and letting $dt \to 0$,

a constant, complex for complex $z$. Evaluating $a(t)$ at the 64 points $z = e^{2\pi i k/64}$ on the unit circle and applying the discrete Fourier transform returns $\Pr(N_t \equiv k \bmod 64)$ for every $k$ at once.

These equal the probabilities $\Pr(N_t = k)$ up to the chance of 64 or more events. Given the regime path the count is Poisson with mean at most 8, so that chance is below $2\times10^{-35}$.

The expansion

Let the regime switch at rate $\lambda$ in each direction, and write

Here

The average $\bar g$ is what the averaged model sees. The half-difference $\tilde g$ measures how much the two regimes disagree.

First order

To first order in $\varepsilon$, at fixed $t > 0$,

with the upper sign for a start in regime 1. If $\tilde g(0)=0$, the remainder is also uniform for $0\leq t\leq T$. Otherwise the uniform composite formula replaces $\tilde g(t)$ in the bracket by

and includes $-\varepsilon\tilde g(0)\int_0^t\tilde g(s)e^{-2s/\varepsilon}ds/2$ in the exponent. Its remainder is $O(\varepsilon^2)$ uniformly on every fixed bounded maturity interval. The proof and the exact crossover scale are on the initial-layers page. The three outer factors have separate meanings, explained on the idea page.

• The factor below is the averaged model: $$e^{\int_0^t \bar g} .$$
• The term below is half the second cumulant of the fluctuating integral of $g$: $$\frac\varepsilon2\int_0^t\tilde g^{\,2} .$$ The square carries no conjugate, and here $g$ is complex whenever $z$ is, and real on the real axis, so the term is complex and is not a variance. For real $g$ it is half the variance, as on the idea page. It is the same for both starting regimes.
• The bracket is the memory of the starting regime. The switch stays near its starting state for about $1/(2\lambda)$, which shifts the integral by $\pm\tilde g/(2\lambda)$.

Every order

Add and subtract the two equations for $a_1$ and $a_2$. The mean $m = \tfrac12(a_1 + a_2)$ and the normalized difference $\omega = (a_1 - a_2)/(a_1 + a_2)$ then satisfy

The first equation has a large restoring term $-2\lambda\omega$: the switch pulls the two regimes together at rate $2\lambda$, while $\tilde g$ pushes them apart. So $\omega$ stays close to the value that balances the two, and expanding in $\varepsilon$ finds that value order by order. Write

and match powers of $\varepsilon$:

This is the two-state case of the solvability argument: at each order the restoring term can be inverted, which here is a division by $2\lambda$.

The series does not satisfy $\omega(0) = 0$.

An initial layer of width about $1/(2\lambda)$ repairs the start and is solved order by order in the same way. Its leading term adds $\mp\tfrac12\varepsilon\,\tilde g(0)\,e^{-2\lambda t}$ to the first-order bracket, which returns the bracket to $1$ at $t = 0$; when $\tilde g(0) \ne 0$ the first-order formula above is therefore accurate to $O(\varepsilon^2)$ only at fixed $t > 0$. The engine page gives the version for any number of regimes.

For the count, the shared first-order term multiplies the generating function by

which adds $t(\ell_1 - \ell_2)^2/(4\lambda)$ to the variance.

The generating function in closed form

The forcing $g_i = (z - 1)\ell_i$ is constant, so the generating function is exact:

with

Differentiating the expansion at $z = 1$ gives the mean and variance of the count. With $L = 1 - e^{-2\lambda T}$ and the upper sign for a start in regime 1,

These two formulas have no higher-order terms. Given the regime path the count is Poisson with mean $\Lambda = \int_0^T \ell_{y_s}\,ds$, so

and integrating the covariance of the two-state regime, which is $e^{-2\lambda|t - s|} - e^{-2\lambda(t + s)}$ in the $\pm1$ coding, gives both lines exactly.

The variance exceeds the mean by

At fixed $T > 0$ the first term leads. It is the Green–Kubo overdispersion. The second term is the memory of the starting regime, and at horizons short against $1/\lambda$ it nearly cancels the first.

At $\ell = (8, 1)$, $T = 1$, $\lambda = 10$ and a start in the busy regime the excess is 1.13313, against 1.22500 from the first term alone:

The numerical column evaluates the generating function at the 64 roots of unity, as described above, and the distribution itself follows the same way. A Monte Carlo check over 400,000 regime paths, sampled exactly from exponential holding times, gives each path a Poisson law for the count.

It agrees with the moments and with every probability up to $k = 15$ within two standard errors. Certificate: verify_option_mc.py.

Factorial cumulants remove Poisson noise

The preceding identity extends to every order and any stochastic intensity whose transform is finite near the origin. Define the factorial cumulant generating function by

Conditional Poisson sampling gives the exact identity

Therefore the $r$th factorial cumulant of $N_T$ is the ordinary $r$th cumulant of $\Lambda_T$. In particular,

Beyond the mean, these combinations strip out Poisson shot noise and leave only fluctuations of the integrated rate. For several counts that are conditionally independent given cumulative intensities $\Lambda_1,\ldots,\Lambda_d$, the joint identity is

Thus mixed factorial cumulants identify the joint cumulants of the integrated rates exactly. For two streams, writing $\kappa_{rs}$ for the ordinary joint cumulant with $r$ copies of $N_1$ and $s$ copies of $N_2$, the signed-Stirling conversion through bidegree $(2,2)$ is

The transform is applied separately in each coordinate. Conditional independence is essential; without it, shared-event shot noise is not removed.

What the count law identifies

At the level of the entire distribution there is a stronger, moment-free result. For a Cox count with nonnegative integrated intensity $\Lambda$,

Hence the full count law determines the Laplace transform of $\Lambda$ on $[0,1]$. Two count laws can agree only if the corresponding Laplace transforms agree on $(0,1)$; analyticity on the positive half-plane and uniqueness of Laplace transforms then imply equality of the mixing laws.

Thus the mixed-Poisson family is identifiable, with no moment-determinacy assumption. This is the Poisson-mixture case of Teicher's identifiability theorem.

What it does not identify: the arrow of latent time

Identifying the law of $\Lambda_T$ is not the same as identifying the hidden generator. Let $Y$ be an irreducible finite-state chain with row generator $Q$, stationary law $\pi$, and $\Pi=\operatorname{diag}(\pi)$. Its stationary time reversal has generator

For state rates $\ell_i\ge0$, put $D=\operatorname{diag}(\ell)$ and $\Lambda_T=\int_0^T\ell_{Y_s}\,ds$. Reversing a path preserves every occupation time, so under a stationary start $Q$ and $Q^\leftarrow$ give exactly the same law of $\Lambda_T$ and hence exactly the same Cox count law for every $T$. This is a genuine non-identification whenever $Q\ne Q^\leftarrow$; the underlying time-reversal construction is classical, going back to Nagasawa (1964).

The endpoint-resolved Feynman–Kac matrices make both the theorem and its boundary explicit. Define

Because $D$ commutes with $\Pi$,

Multiplication on the left by $\pi$ and on the right by $\mathbf1$ erases this weighted transpose, proving

The same proof covers several conditionally independent count streams by replacing $(z-1)D$ with a sum of diagonal state-rate potentials. It also covers state-dependent compound-Poisson marks by using their diagonal probability-generating exponents. Thus neither additional terminal count coordinates nor arbitrary contemporaneous state marks recover the arrow of time if all observations remain occupation based.

The stationarity and aggregation qualifications are essential. With a known start $i$, the transform is $e_i^\top K_Q\mathbf1$, which generally differs after reversal.

Recording the terminal state also retains the weighted-transpose relation above, and ordered event-time or path data can retain temporal direction. The theorem concerns terminal occupation-based counts, not those richer experiments.

Ordered windows recover the arrow of time

A minimal enrichment already breaks that non-identification. Let $N_a^-$ count stream $a$ on $[0,A]$ and $N_b^+$ count stream $b$ on $[A+\Delta,A+\Delta+B]$, where stream $a$ has state rate $\ell^{(a)}$.

Write $F_{ia}=\ell_i^{(a)}-\pi\ell^{(a)}$, $P_t=e^{tQ}$ and $H_T(Q)=\int_0^T P_s\,ds$. Conditional Poisson increments on disjoint windows are independent, so there is no shot-noise term across the boundary and the entire early/late covariance matrix is exactly

Since $e^{tQ^\leftarrow}=\Pi^{-1}e^{tQ^\top}\Pi$, time reversal transposes this matrix:

Thus the antisymmetric part is an exact arrow-of-time statistic. It vanishes for reversible chains, but it need not vanish for irreversible ones.

This conclusion uses ordered, stream-labelled windows; summing the windows returns to the occupation-only experiment and loses the direction. This is a finite-window, time-domain counterpart of the second-order point-process analysis initiated by Bartlett (1963).

The fast-switching scale is a boundary effect. Put $Q_m=mQ_0$, take adjacent fixed windows $A,B>0$, and let $R=-Q_0^\#$ be the zero-mean potential operator. Then

More generally, for a microscopic gap $\Delta=\delta/m$ the middle factor survives and the limit is $F^\top\Pi R e^{\delta Q_0}RF$; for a fixed positive gap the covariance instead decays exponentially. This $m^{-2}$ ordered-boundary signal is distinct from the $m^{-1}$ Green–Kubo variance accumulated inside one macroscopic window.

For the four-state irreversible certificate with $A=0.7$ and $B=1.1$, the off-diagonal antisymmetric signal is 0.022655140619. Reversal agrees with the transpose to 3.1e-17; an independent mixed derivative of the sequential count PGF agrees with the semigroup covariance to 2.2e-16. At $m=40$, the scaled covariance is within 9.4e-16 of $F^\top\Pi R^2F$.

The certificate uses a nonreversible four-state generator whose maximum entrywise difference from its reverse is 1.875637241749. Across four PGF arguments, stationary transforms agree within 4.4e-16 and the weighted endpoint identity within 3.9e-16.

An independent count-resolved master equation agrees within 2.8e-16, with omitted mass below 3.1e-64. In contrast, a known start in state 1 gives count-law total variation 0.068465756461.

Which count laws are mixed Poisson?

Injectivity answers uniqueness after a mixed-Poisson representation is known to exist. Existence itself has an exact test.

Let $p_n=\Pr\{N=n\}$ and $q_n=n!p_n$. Then $N$ is mixed Poisson if and only if, for every $r\ge0$, both Hankel matrices

are positive semidefinite. This is precisely the classical Stieltjes moment criterion; its connection with mixed Poisson laws and factorial moments is also discussed by Kuba and Panholzer.

The proof is short and also explains the factorial weighting. If $p_n=\int e^{-\lambda}\lambda^n/n!\,\mu(d\lambda)$, then

so both matrix families are Gram matrices of polynomials in $L^2(\rho)$ and hence are positive semidefinite. Conversely, the two Hankel conditions give a measure $\rho$ on $[0,\infty)$ with moments $q_n$. Since $\sum_nq_n/n!=\sum_np_n=1$, monotone convergence gives $\int e^\lambda\rho(d\lambda)=1$.

Therefore $\mu(d\lambda)=e^\lambda\rho(d\lambda)$ is a probability measure and reproduces $p_n$. The same exponential moment makes $\rho$ moment-determinate, so the mixing law is unique.

This yields a nested falsification hierarchy. The first two nontrivial minors require

but finitely many such inequalities are not sufficient. For example, let $B$ have probabilities $(0.3,0.1,0.1,0.2,0.3)$ on $\{0,1,2,3,4\}$, let $\Pr\{G=n\}=2^{-n-1}$, and take $N$ to have law $0.9\mathcal L(B)+0.1\mathcal L(G)$. This law has full support and is overdispersed:

Its two order-one Hankel determinants are 0.0523750 and 0.0864875, yet

which proves it is not mixed Poisson. As a positive control, the certificate constructs a four-atom mixing law; after diagonal normalization, the smallest eigenvalues of $H_2^{(0)}$ and $H_2^{(1)}$ are 0.03337128 and 0.02082002. Thus overdispersion is only the first coarse necessary condition, while the full Hankel hierarchy is necessary and sufficient.

The multivariate statement is identical under conditional independence:

Finite cumulant lists are fundamentally weaker. Fix any order $k$, set $n=k+1$ and $x_j=j+1$ for $j=0,\ldots,n$. For any $0<\delta<2^{-n}$, define two strictly positive probability laws by

The $n$th finite-difference identity gives

Consequently the two intensities have the same first $k$ moments and cumulants, so their mixed-Poisson counts have the same first $k$ factorial cumulants. They are nevertheless different count laws because

For the certificate's $k=4$ construction, the maximum discrepancy among the first four factorial cumulants is 0.00e+00, while the zero-count probability differs by 0.000928207565 and the numerically summed total-variation distance is 0.003607864895, with omitted contribution below 6.68e-60. Full-law identification, all-order analytic identification, and finite-order cumulant identification are therefore distinct claims.

The sharp finite-atomic exception

A sparsity assumption changes the last conclusion completely. Suppose the mixing law has at most $r$ atoms,

Then factorial cumulants of orders $1,\ldots,2r-1$ identify $\mu$, and this order is sharp without additional separation or weight assumptions. Indeed, factorial cumulants of the count are ordinary cumulants of $\Lambda$, and the triangular moment-cumulant relations recover $m_k=\mathbb E[\Lambda^k]$ for $0\le k\le2r-1$, with $m_0=1$.

To reconstruct the law, first read $s$ as the rank of $H_{r-1}=(m_{i+j})_{i,j=0}^{r-1}$. For the actual support size,

Let $P(x)=x^s+c_{s-1}x^{s-1}+\cdots+c_0$ be the monic support polynomial. Since $P(x_j)=0$, its coefficients are the unique solution of the Hankel system

The roots of $P$ are the atoms $x_j$, after which the first $s$ moment equations form a nonsingular Vandermonde system for the weights. This is the classical annihilating-polynomial mechanism behind Prony's method.

The count $2r-1$ cannot be reduced. On the $2r$ nodes $x_j=j+1$, split the signed coefficients

into their positive and negative parts and normalize each by $2^{2r-2}$. Each part is an $r$-atomic probability law. The finite-difference identity makes their moments, and hence cumulants, identical through order $2r-2$, while their order-$(2r-1)$ moment gap is

For $r=4$, the certificate recovers a four-atom law with maximum support and weight errors 1.76e-12 and 5.96e-13. Its exact sharpness pair agrees through order six and has order-seven moment and cumulant gap -78.75. This result is purely about exact identification: Hankel and Vandermonde systems can become arbitrarily ill-conditioned as atoms collide or weights vanish, so it does not assert uniform stable recovery.

Exact local conditioning geometry

The instability boundary has an exact algebraic description. For an exactly $r$-atomic law, eliminate the last weight by $w_r=1-\sum_{j<r}w_j$ and define the square parameter-to-moment map

On the ordered interior $w_j>0$ and $x_1<\cdots<x_r$, its Jacobian satisfies

The moment-to-cumulant transformation is unit triangular, so the parameter-to-factorial-cumulant map has the same determinant.

To prove the identity, temporarily retain all $r$ weights and include $m_0$. The $j$th weight and atom columns of the full Jacobian are respectively

After interlacing these columns, this is a weighted confluent Vandermonde matrix, whose determinant is the fourth power of the ordinary Vandermonde product times $\prod_jw_j$. Changing coordinates from the weights to $(w_1,\ldots,w_{r-1},m_0)$ and expanding along $m_0=1$ leaves the displayed reduced determinant. The classical inverse and conditioning analysis of these matrices goes back to Gautschi (1962).

The inverse-function theorem now gives a locally analytic inverse everywhere in the ordered interior. It also yields an explicit differential bound.

Put $d=2r-1$, suppose $0\le x_j\le R$, $w_j\ge w_*$, and $x_{j+1}-x_j\ge\delta$, and set $A_R=\max\{1,R^{2r-1}\}$. Every entry of $D\Phi$ is bounded by $dA_R$, so

Thus separated atoms with weights bounded away from zero have uniformly bounded infinitesimal inverse sensitivity. The determinant formula also pinpoints both degeneracies: it vanishes linearly with a disappearing weight and to fourth order for each colliding pair. This is a local parameter-stability statement, not a claim of total-variation stability for unrestricted mixing measures.

The certificate evaluates the reduced Jacobian directly for two through five atoms. Its maximum relative determinant error is 2.03e-11; a centered directional derivative agrees within 8.83e-09; and after removing the noncolliding factors, a shrinking pair has measured volume-collapse order 4.000000001.

Identification is not stable inversion

Injectivity is qualitative. It does not make recovery of an unrestricted mixing law stable in total variation. For two point-mass mixing laws $\delta_a$ and $\delta_b$,

whereas their count laws are $\operatorname{Pois}(a)$ and $\operatorname{Pois}(b)$. Their Hellinger affinity is

so

For $0<a<b$, the likelihood ratio $p_a(k)/p_b(k)=e^{b-a}(a/b)^k$ crosses one once. Therefore the exact total variation is

Thus even on a fixed compact intensity interval, taking $b\to a$ leaves mixing-law total variation equal to one while count-law total variation tends to zero. More precisely, at a noninteger $a$ and for sufficiently small $h>0$, the likelihood-ratio crossing is $k=\lfloor a\rfloor$, and

There is no global Wasserstein rescue without an intensity bound. For integer $m$, the crossing of $\operatorname{Pois}(m)$ and $\operatorname{Pois}(m+1)$ is exactly $k=m$, so

although $W_1(\delta_m,\delta_{m+1})=1$. At $a=4.5$ and $h=10^{-5}$ the certificate gives $d_{\rm TV}/h=$ 0.189807515050, against the limit 0.189807620540. At $m=10^6$ it gives $\sqrt m\,d_{\rm TV}=$ 0.398942180666, against $1/\sqrt{2\pi}=$ 0.398942280401.

This does not contradict identifiability. It says that finite-sample recovery needs a weaker loss and/or structural restrictions such as bounded support and smoothness.

The nonparametric Poisson-mixture estimation theory of Roueff and Rydén (2005) makes such regularity assumptions explicit. For a fixed finite-state intensity model, $\Lambda_T$ is bounded; that removes the escaping-mass example but not the total-variation discontinuity created by moving atoms.

A finite-sample impossibility theorem

The discontinuity implies more than the absence of a convenient inverse bound. Suppose $N_1,\ldots,N_n$ are iid mixed-Poisson observations with unknown mixing law $\mu$ supported on a fixed compact interval containing more than one point. For any estimator $\widehat\mu_n$, including a randomized estimator,

Thus the unrestricted mixing law is not uniformly consistently estimable in total variation, even though it is identified by the population count law and even though its support is compact.

The proof is a two-point coupling argument. For any two mixing laws $\mu_0,\mu_1$, let $P_i$ be the corresponding one-count laws and put $R_i=\mathbb E_i d_{\rm TV}(\widehat\mu_n,\mu_i)$. A maximal coupling of $P_0^{\otimes n}$ and $P_1^{\otimes n}$, using the same estimator randomization when the samples agree, and the triangle inequality give

Choose distinct point masses $\mu_0=\delta_a$ and $\mu_1=\delta_b$ inside the interval. Their distance is one, while the product count distance tends to zero as $b\to a$. Letting $b$ approach $a$ proves the displayed lower bound for every fixed $n$.

There is also an exact quantitative certificate. Take $a=4.5$, $b_n=a+1/n$, and even $n$.

Under the two hypotheses, the sufficient statistic $\sum_iN_i$ has laws $\operatorname{Pois}(na)$ and $\operatorname{Pois}(na+1)$; conditional on the sum, the allocation among the $n$ observations is the same multinomial law. Hence

At $n=10^6$, $\sqrt n$ times the product distance is $0.188063184068$, versus $1/\sqrt{2\pi a}=0.188063194516$. The lower bound concerns total-variation recovery of an unrestricted measure. It does not rule out weak-loss consistency, parametric finite-state recovery under separation, or density estimation on smoothness classes; those are precisely the kinds of restrictions used by Roueff and Rydén.

A local Wasserstein lower bound

Compact support does make Wasserstein loss qualitatively weaker than total variation, but it does not permit a uniformly faster-than-parametric rate. Let $\mathcal M_{[4,5]}$ be all probability laws on $[4,5]$ and define

For any $a$ in the interior of the interval and any fixed $h>0$, compare $\mu_{0,n}=\delta_a$ with $\mu_{1,n}=\delta_{a+h/\sqrt n}$. They lie in the class for all sufficiently large $n$ and have $W_1(\mu_{0,n},\mu_{1,n})=h/\sqrt n$. The same two-point coupling inequality used above gives

This is a restriction to the point-mass submodel, so it is a valid lower bound for the full class. It is the classical two-point testing reduction of Le Cam, here with the distance and testing affinity evaluated exactly.

Put $z=h/(2\sqrt a)$. The one-crossing formula has cutoff

Applying the normal limit to the two Poisson distribution functions at this cutoff yields

and therefore

The right side is optimized by $h=2\sqrt a\,z_*$, where $z_*$ is the unique positive root of $\Phi(-z)=z\phi(z)$. At $a=4.5$,

The final column tends to $0.721126760493$. This proves only a lower bound: it rules out uniform $o(n^{-1/2})$ Wasserstein-1 recovery, even on the point-mass submodel, but it does not claim that the unrestricted mixing class has an $O(n^{-1/2})$ estimator.

A nonparametric logarithmic obstruction

The unrestricted compact class is substantially harder than its point-mass submodel. More generally, for $a\ge0$ and $B>0$, put $\mathcal M_{[a,a+B]}$ for all laws on that intensity interval, and let $R_n^{(1)}(a,B)$ denote the same minimax risk with this class. Then

Thus no estimator has worst-case $W_1$ risk $o((\log\log n)/\log n)$ on the unrestricted class.

The construction is explicit. Fix an integer $L\ge2$, let $J\sim\operatorname{Bin}(L,1/2)$, put $h=B/L$, and define

The finite-difference identity

shows that the two laws have identical first $L-1$ moments. Yet on $[a+kh,a+(k+1)h)$ their CDF difference is

so summing the absolute areas gives the exact separation

The induced mixed-Poisson count laws are nevertheless exponentially close in $L$. To see this without an analytic inversion bound, write a count as

where $Z\sim\operatorname{Pois}(a)$, the $B_i$ are fair Bernoulli variables conditioned on even or odd parity, and the independent $X_i\sim\operatorname{Pois}(h)$. Whenever some $X_i=0$, flipping the first corresponding $B_i$ changes parity without changing $N$ and maps the uniform even-parity law bijectively to the uniform odd-parity law. Therefore, for the one-count laws $P_{L,+},P_{L,-}$,

and product coupling gives

Le Cam's two-point metric inequality now yields, for every $L\ge2$,

For fixed $\epsilon>0$, take $L=\lceil(1+\epsilon)\log n/\log\log n\rceil$. Since $n(1-e^{-B/L})^L\le n(B/L)^L\to0$, the scaled lower limit is at least $B/[2(1+\epsilon)]$; letting $\epsilon$ decrease to zero proves the boxed result.

The certificate below first checks exact moment matching, exact $W_1$, and the count-law coupling for $[a,a+B]=[4,5]$. It then optimizes the displayed finite-$n$ lower bound over integer $L$.

The scaled certificate approaches its proved asymptotic lower constant $1/2$ slowly. This ordinary, unsmoothed $W_1$ obstruction is consistent with the smoothness-class theory of Roueff and Rydén. For support $[0,B]$, Miao et al. prove the matching upper order for the NPMLE and a matching minimax lower order; Lim and Han obtain a nearly root-$n$ rate only after replacing ordinary transport by Gaussian-smoothed optimal transport.

A matching moment-estimator upper rate

The upper order can also be recovered by a direct estimator on any fixed interval $[a,a+B]$. Put $M=a+B$ and, from iid counts $N_1,\ldots,N_n$, define the unbiased normalized factorial-moment estimates

Choose any probability law on $[a,M]$ minimizing the largest moment residual,

A minimizer exists by weak compactness; the compact-valued argmin correspondence has a Borel graph, so fix any measurable selection. If $\epsilon_L=\max_{k\le L}|\widehat m_k-m_k|$, comparison with the true law gives a discrepancy of at most $2\epsilon_L$ in every fitted moment.

For every $n,L\ge1$, the following nonasymptotic bound is explicit up to the universal Jackson constant $C_J$:

To prove it, use Kantorovich duality and subtract the value of each Lipschitz test function at $a$. Extend it constantly to $[0,a]$, rescale to $[0,1]$, and apply Jackson's polynomial approximation theorem. A degree-$L$ polynomial approximates the test function within $C_JM/L$.

In the shifted Chebyshev basis $q_j(x)=T_j(2x-1)$, the recurrence $q_{j+1}=(4x-2)q_j-q_{j-1}$ implies that the monomial coefficient $\ell^1$ norm of $q_j$ is at most $7^j$. The Chebyshev coefficients of a bounded polynomial are at most twice its sup norm, so the approximant's monomial coefficient norm is at most $3(B+C_JM)7^L$.

The stochastic term follows from the exact falling-factorial product identity

Since $N\mid\theta\sim\operatorname{Pois}(\theta)$ and $\theta\le M$, it gives

Combining the approximation and moment errors proves the displayed bound. For any fixed $0<c<1$, take $L_n=\lfloor c\log n/\log\log n\rfloor$.

The logarithm of the stochastic factor is $-\tfrac12(1-c)\log n+o(\log n)$, while the approximation term is $O(\log\log n/\log n)$. Consequently,

Together with the preceding lower theorem, this proves that the unrestricted compact minimax rate is exactly $\Theta((\log\log n)/\log n)$ for every fixed $a\ge0$ and $B>0$. The certificate checks the factorial second-moment identity through order eight to relative error 2.89e-15 and verifies the coefficient recurrence. For $M=5$ and $c=1/2$, the derived stochastic factor decays as follows.

The compact inverse has an optimal logarithmic modulus

The statistical upper bound also yields a population stability theorem. For a mixing law $\mu$ on $[a,a+B]$, write $P_\mu$ for its one-count mixed-Poisson law and define

Fix one measurable version of the moment estimator above. Maximally couple $n$ iid counts from $P_\mu$ and $P_\nu$, and apply the same estimator to the coupled samples. They agree with probability at least $(1-\delta)^n$.

On the complementary event, two estimator outputs supported on $[a,a+B]$ are at Wasserstein distance at most $B$. The triangle inequality therefore gives

where $R_{n,L}$ is the preceding worst-case estimator bound. In particular, with $M=a+B$ and $A_M=\max(1,M^{-1})$,

This holds for every positive integer $n,L$. Taking $n=\lfloor\delta^{-1/2}\rfloor$ and $L=\lfloor c\log n/\log\log n\rfloor$ for any fixed $0<c<1$ makes the last two terms negligible relative to $1/L$, and proves

The order cannot be improved. For the even/odd binomial pair used in the minimax lower bound,

Hence $\omega_{a,B}(\delta_L)\ge B/L$. More generally, for every sufficiently small $\delta$, choosing $L$ of order $\log(1/\delta)/\log\log(1/\delta)$ makes $\delta_L\le\delta$ and gives the reverse order. Therefore

This is a positive compact-support stability statement in $W_1$, not in total variation. It complements the minimax theorem of Miao et al.: the same logarithmic ill-posedness governs both population inversion and unrestricted statistical recovery.

On $[4,5]$, the certificate's lower scaled separation reaches 0.781716448592 at $L=1024$ on its slow approach to one. At $\log(1/\delta)=10{,}000$, the upper construction selects $L=293$; the logarithms of its stochastic and sample-disagreement terms are -886.085546 and -5000.000000.

The point-mass submodel has the parametric rate

The preceding lower bound is rate-sharp on the submodel that generated it. Define

When the truth is $\delta_\lambda$, the observations are iid $\operatorname{Pois}(\lambda)$. Put $S_n=\sum_iN_i$, project $S_n/n$ onto $[4,5]$, and put a point mass at the projected value. Since projection cannot increase distance to $\lambda$ and $S_n\sim\operatorname{Pois}(n\lambda)$,

The equality follows by splitting the centered Poisson variable at its mean and using $k p_\mu(k)=\mu p_\mu(k-1)$; the inequality is Cauchy–Schwarz. Combining this estimator with the local two-point lower bound gives the rigorous sandwich

So $R_{n,\delta}^{(1)}=\Theta(n^{-1/2})$. At the certificate point $\lambda=4.5$, projection changes the interior risk only by exponentially small tails, and the unprojected sample-mean benchmark has the sharper pointwise limit

This sample-mean upper bound applies only to the one-parameter family of point masses. An arbitrary law in $\mathcal M_{[4,5]}$ cannot be estimated by reducing it to its mean; the moment estimator above attains the slower unrestricted logarithmic rate.

What common shocks add

The correction is exact. Index independent Poisson event streams by the nonempty subsets $A\subseteq\{1,\ldots,d\}$ of coordinates that each event increments.

Conditional on their cumulative intensities $\Lambda_A$, put $N_j=\sum_{A\ni j}C_A$. Then

For two counts, write $\Lambda_0$ for the common-event intensity and $A=\Lambda_1+\Lambda_0$, $B=\Lambda_2+\Lambda_0$, $C=\Lambda_0$. Expanding the preceding identity through arbitrary bidegree $(r,s)$ gives the all-order identity

Here the brackets denote repeated arguments of the joint cumulant; when only one argument remains, it is its expectation. The coefficient counts which $k$ derivatives in each coordinate strike the bilinear common-event term. In particular, through bidegree $(2,2)$,

Thus mixed factorial cumulants identify marginal-intensity dependence only after the common-shock terms are modeled or ruled out. The sharp counterexample has deterministic $C=cT$ and no idiosyncratic intensity: $N_1=N_2\sim\operatorname{Poisson}(cT)$.

The cumulative intensities have zero covariance, but $\kappa^{(F)}_{11}=cT$. At $c=0.8$ and $T=1.3$, the certificate obtains $1.04000000$.

Arbitrary integer marks

The unit-jump restriction is unnecessary. Let each common event carry an iid nonnegative integer mark $J=(J_1,J_2)$, independent of the cumulative intensities and of the other marks.

Conditional on the common cumulative intensity $C$, the number of marked events is Poisson with mean $C$. If $\Lambda_1,\Lambda_2$ are the idiosyncratic cumulative intensities, then the exact factorial cumulant generating function is

Put $\mu_{ab}=\mathbb E[(J_1)_a(J_2)_b]$, where $(j)_a$ is a falling factorial, and define the integrated-rate variables

Take $r$ labelled symbols of type 1 and $s$ of type 2. For a block $B$ of a set partition, let $a_B,b_B$ be its numbers of the two symbol types. The all-order marked theorem is

This is the multivariate logarithmic Faà di Bruno formula applied to $H$. It requires the displayed falling-factorial mark moments and the corresponding integrated-rate cumulants to be finite. For example,

When $J=(1,1)$ almost surely, every $\mu_{ab}$ vanishes except those with $a,b\le1$. Only singleton blocks and disjoint type-1/type-2 pairs survive, and the partition formula reduces exactly to the preceding ${r\choose k}{s\choose k}k!$ matching formula. Thus larger or asymmetric common marks do not invalidate factorial-cumulant identification; they change the shared-shot-noise correction through their factorial moments.

Regime-dependent mark laws

The iid mark law may be replaced by a law that depends on the hidden regime at the event time. Write $c_i$ for the common-event rate and $\mu_{ab}(i)=\mathbb E[(J_1)_a(J_2)_b\mid Y_t=i]$. Conditional on the complete regime path, the exact factorial cumulant generating function becomes

where $M_i(t,u)=\mathbb E[(1+t)^{J_1}(1+u)^{J_2}\mid Y=i]$. Therefore the same partition theorem holds after replacing the variables above by

This is a finite-state Markov additive process in the sense of Çinlar (1972). The condition is pathwise: marks are conditionally independent given the regime path and their law at an event depends only on the contemporaneous regime. It does not cover dependence between different event marks beyond that hidden-state dependence.

The certificate obtains the factorial-count and integrated-intensity columns independently from the regime/count master equation and a polynomial Feynman–Kac hierarchy. It checks the univariate identity through order four and the bivariate identity through bidegree $(2,2)$ for a nonreversible three-state chain. The mixed calculation agrees within $5.3\times10^{-14}$, with omitted count mass bounded by $1.6\times10^{-47}$.

A separate master equation with simultaneous $(1,1)$ jumps checks the common-shock formula through every mixed order $(r,s)$ with $1\leq r,s\leq3$ within $7.9\times10^{-12}$; at $(3,3)$ the two independent values are 3.05155621 and 3.05155621. A second marked master equation uses six mark vectors, including $(3,1)$ and $(1,3)$, and agrees with the marked partition theorem through bidegree $(3,3)$ within $9.6\times10^{-10}$; the independent values of $\kappa^{(F)}_{33}$ are 1.41746379 and 1.41746379, with omitted mass below $8.7\times10^{-49}$. Substituting the degenerate mark $(1,1)$ recovers the matching formula within $1.1\times10^{-11}$.

A third master equation uses a different six-point mark distribution in every regime. It agrees with the state-dependent partition formula through bidegree $(3,3)$ within $6.5\times10^{-10}$; at $(3,3)$ the independent values are -1.36187430 and -1.36187430, and omitted mass is below $1.4\times10^{-44}$. The partition calculation follows the joint-cumulant method of Leonov and Shiryaev (1959).

The certificate also measures second-order decay of the univariate omitted boundary term. The conditioning argument is the defining construction of a Cox process, specialized here to the Markov-modulated Poisson process. The common-component construction is the classical bivariate Poisson model of Holgate (1964); the arbitrary-mark construction is a multivariate compound-Poisson law in the sense of Feller (1943).

What the count identifies

Given the regime path the count is Poisson with mean $\Lambda_T = \int_0^T \ell_{y_s}\,ds$, so $\mathbb E[(1 + t)^{N_T}] = \mathbb E[e^{t\Lambda_T}]$: the factorial cumulant generating function of the count is the cumulant generating function of the integrated rate. Every factorial cumulant of the count equals the corresponding cumulant of $\Lambda_T$:

These combinations strip out the Poisson shot noise and keep the fluctuations of the rate. The first is the overdispersion above, which is $2TK$ with $K$ the Green–Kubo integral of the rate up to a remainder bounded uniformly in $T$, as on the cumulants page; for two regimes the remainder is at most $\varepsilon^2\tilde\ell^{\,2}/2$, the size of the memory term.

The third and fourth are new: they measure the skewness and kurtosis of the integrated rate, which a test for switching in count data would examine, as the likelihood page examines the residuals of a rate. For a stationary start with two symmetric regimes the third vanishes, and the fourth is negative because the rate is bounded. At $T = 1$, computing the factorial cumulants from the count law at the roots of unity and the cumulants of $\Lambda_T$ from its polynomial moment system, the two agree to the rounding of the count law:

Certificate: verify_counts.py.

Results

Moments of $N_1$, starting in the busy regime, from the numerical solution:

The largest error over all probabilities after each order: