ACID: A Low Dimensional Characterization of Markov-Modulated and Self-Exciting Counting Processes

Mark Sinzger-D'Angelo and
Heinz Koeppl

This paper was presented in part at the 2020 IEEE International Symposium on Information Theory, virtual, Los Angeles, California, USA.
The work was supported by the European Research Council (ERC) within the Consolidator Grant CONSYN (grant agreement no. 773196).
M. Sinzger-D'Angelo and H. Koeppl are with the Department of Electrical Engineering and Information Technology, Technische Universität Darmstadt (email: {mark.sinzger, heinz.koeppl}@tu-darmstadt.de). H. Koeppl is with the Centre for Synthetic Biology, Technische Universität Darmstadt.

The conditional intensity (CI) of a counting process $\out_\ti$ is based on the minimal knowledge $\sigmalg^\out_\ti$, i.e., on the observation of $\out_\ti$ alone. Prominently, the mutual information rate of a signal and its Poisson channel output is a difference functional between the CI and the intensity that has full knowledge about the input. While the CI of Markov-modulated Poisson processes evolves according to Snyder's filter, self-exciting processes, e.g., Hawkes processes, specify the CI via the history of $\out_\ti$. The emergence of the CI as a self-contained stochastic process prompts us to bring its statistical ensemble into focus. We investigate the asymptotic conditional intensity distribution (ACID) and emphasize its rich information content. We assume the case in which the CI is determined from a sufficient statistic that progresses as a Markov process. We present a simulation-free method to compute the ACID when the dimension of the sufficient statistic is low. The method is made possible by introducing a backward recurrence time parametrization, which has the advantage to align all probability inflow in a boundary condition for the master equation. Case studies illustrate the usage of ACID for three primary examples: 1) the Poisson channels with binary Markovian input (as an example of a Markov-modulated Poisson process), 2) the standard Hawkes process with exponential kernel (as an example of a self-exciting counting process) and 3) the Gamma filter (as an example of an approximate filter to a Markov-modulated Poisson process).

Conditional intensity, Poisson channel, mutual information rate, average sojourn time constraint, Snyder filter, backward recurrence time, marginal simulation, Hawkes process.

§ INTRODUCTION

The study of the Poisson channel dates back to the 70s. While originally introduced as a model of optical communication [1], its setting is ubiquitous in various fields of research for modeling event counts: An input signal $\msg_\ti$, possibly corrupted by dark current, is partially observed via point observations. In the terminology of Cox [2] and Snyder [3], the channel output $\out_\ti$ is a doubly stochastic Poisson process with feedback-less intensity $\intensity_\ti = \intensityFct(\msg_\ti)$.

Much is already known about the Poisson channel. In particular, its capacity under amplitude constraint $\msg_\ti \in [0,\energy]$ was determined by Kabanov [4] shortly after its introduction. However, the capacity-achieving input distribution is found to switch infinitely fast between $0$ and $\energy$ [4, 5]. This flaw in terms of physical interpretability was addressed by adding constraints on the bandwidth [6, 7, 8]. While reducing the class of input processes this necessitates to review the task of computing the MI.

Besides taking its definition

\begin{equation}
\I(\msg_{[0,T]}, \out_{[0,T]}) := \E\left[ \ln \frac{\dd \mu^{\msg,\out}_T}{\dd \mu^\msg_T \times \mu^\out_T} \right]
\label{def:MI}
\end{equation}

as a departure for the computation, an expression found to be useful was obtained by Liptser [9]

\begin{equation}
\I(\msg_{[0,T]}, \out_{[0,T]}) = \int_0^T \E[\phi(\intensity_\ti) - \phi(\hat \intensity_\ti)] \dd \ti,
\label{eq:Liptser_time}
\end{equation}

where $\phi(z) = z \ln z$ and $\hat \intensity_\ti = \E[\intensity_\ti\vert \sigmalg^\out_\ti]$. Although other results on links between the MI and conditional estimation received considerable attention [10, 11], the classical expression eq. (<ref>) remains the standard (besides the definition eq. (<ref>)) when searching for ways to compute the MI [12, 13].

The difficulty in computing the MI via eq. (<ref>) clearly lies in the computation of $\E[\phi(\CIntensity_\ti)]$. Different approaches address this. (i) In case of a stationary input $\intensity_\ti$, one can consider

$$ \intensity_\ti \to \intensity_\infty, \quad \CIntensity_\ti \to \CIntensity_\infty$$

in distribution, still allowing conclusions about the mutual information rate

\begin{equation}
\IR(\msg, \out) : = \lim_{T \to \infty}\frac 1 T \I(\msg_{[0,T]}, \out_{[0,T]}) = \E[\phi(\intensity_\infty)] - \E[\phi(\hat \intensity_\infty)].
\label{eq:Liptser}
\end{equation}

(ii) Properties of the function $\phi$ (Lipschitz for bounded $\intensityFct$ [4] or non-negative third derivative [8]) were used for a reduction to second order moments of $\CIntensity_\infty$. (iii) Kabanov, Davis used martingale theory with its rich second order analysis tools.
(iv) The conditioning on a coarser sigma-field [7] or the use of a suboptimal estimator [8] provided upper bounds. (v) Alternatively, one needs to have the knack for choosing tractable input process classes. For example, piecewise-constant input trajectories, whose amplitudes follow a Markov chain were considered [14].

We took the approach (i) combined with (v), restricting ourselves to continuous-time Markov chains $\msg_\ti$ with low number of input states.
Its appearance on the right-hand side of (<ref>) is our main motivation to look at the distribution of $\CIntensity_\infty$, the asymptotic conditional intensity distribution (ACID). With the spotlight being on $\CIntensity_\ti$, we elaborate on its interpretation as conditional intensity in <ref>. A summary of our general contributions is followed by a confined focus on cell biology for interested readers (<ref>). In <ref> we advocate that the study of ACID can be beneficial beyond the information rate.

§.§ Conditional intensity view replaces estimator perspective

The conditional mean $\CIntensity_\ti = \E[\intensity_\ti \vert \sigmalg^\out_\ti]$ is often regarded as an estimator of $\intensity_\ti$. There is a clear justification for this: For each $\ti \geq 0$

\begin{equation}
\E[\intensity_\ti \vert \sigmalg^\out_\ti] = \arg \min_{Z} \E[(\intensity_\ti - Z)^2]
\end{equation}

among all random variables $Z$ that are $\sigmalg^\out_\ti$-measurable. But this interpretation attributes a notion of deficit to $\CIntensity_\ti$: While $\intensity_\ti$ is the true intensity of $\out_\ti$, we obtain merely its optimal yet inaccurate approximation via $\CIntensity_\ti$. We advocate a different perspective throughout the paper.

The process $\out_\ti$ is a self-exciting process with intensity $\CIntensity_\ti$, which has been known at least since the monograph [3] by Snyder. How can $\out_\ti$ have two intensities? How can $\CIntensity_\ti$ be on a par with $\intensity_\ti$?
This puzzle is solved if one incorporates that the intensity $\intensity_\ti$ of a general counting process $\out_\ti$  has an $\sigmalg_\ti$-dependency. This dependency can be expressed in at least two ways: (i) explicitly but not of operational value for mathematical proofs

$$ \intensity_\ti = \E[\dd \out_\ti\vert \sigmalg_\ti], $$

(ii) implicitly, and beneficial for mathematical rigidity [15] by requesting

$$ \out_\ti - \int_0^\ti \intensity_s \dd s$$

to be an $\sigmalg_\ti$-martingale.
It implies that a general counting process $\out_\ti$ can possess several intensities. For a strictly coarser $\mathcal G_\ti \subseteq \sigmalg_\ti$ and corresponding intensities $\hat \intensity_\ti, \intensity_\ti$, Bremaud's innovation theorem provides
$ \hat \intensity_\ti = \E[\intensity_\ti \vert \mathcal G_\ti] $.
This allows the following interpretation of eq. (<ref>). The information rate is not a difference functional between the true intensity and its deficient estimate, but between two intensities that differ in the state of knowledge: $\sigmalg^{\msg,\out}_\ti = \sigma(\msg_s,\out_s: s \le \ti)$ vs. $\sigmalg^{\out}_\ti = \sigma(\out_s: s \le \ti)$. We follow Bremaud [16] and Daley, Vere-Jones [17] in using the term conditional intensity for the $\sigmalg^\out_\ti$-intensity.
It is not to be confused with its other usage as a special case of the Papangelou intensity.

Far from being an estimator with deficits, the conditional intensity $\CIntensity_\ti$ is canonical one among the $\sigmalg^{\out}_\ti$- and $\sigmalg^{\msg,\out}_\ti$-intensities. The $\sigmalg^{\out}_\ti$-intensity can be defined for any counting process $\out_\ti$, with no need to specify an external process.
And it has a standing as being the intensity with respect to the minimal filtration to which $\out_\ti$ is adapted.

The estimator perspective regards $\E[\intensity_\ti \vert \sigmalg^\out_\ti]$ as a function of the history $\{ \out_s \colon 0 \le s \le \ti \}$. This
complicated functional dependence makes it hard
to understand $\CIntensity_\ti \to \CIntensity_\infty$. In contrast, the conditional intensity view offers a self-contained description of $\CIntensity_\ti$. It emancipates from $\out_\ti$ when both have common jumps. This view guides our understanding of $\CIntensity_\ti \to \CIntensity_\infty$.

We restrict our study of ACID to processes $\out_\ti$ for which $\CIntensity_\ti$ is obtained from a piecewise-deterministic Markov process. Within the scope of this work, we review that Markov-modulated Poisson processes belong to this class [18], as do the Hawkes process [19, 20] and approximate filters, introduced in section <ref>. We exploit that the ACID can then be obtained from generator theory of Markov processes[21]. To this end, we contribute a backward recurrence time parametrization (BReT-P), whose advantage is the alignment of all probability inflow. Since our numerical method works via grid discretization, it is limited in the dimension of the piecewise-deterministic Markov process. But conceptually, BReT-P opens the door for new numerical methods that compute the MI without Monte Carlo simulation.
We illustrate how our method helps with screening for system parameters that optimize the information rate under bandwidth-like constraints. Applications beyond information theory complement the case studies.

§.§ Path mutual information in cell signaling

Our particular interest lies in studying the mutual information between time-varying signals in biological systems. We review the state of the art of information theory in cell biology to locate our contribution to the field. Readers mainly interested in the method can skip this paragraph.

It is widely assumed that the optimization of the MI could be an evolutionary strategy to mitigate noise in the cellular signalling via metabolites and in gene regulatory networks [22]. There is mounting evidence that information is encoded in the temporal profile of biomolecules [23, 24, 25], but the exact mechanisms how temporal features result in decision making are unclear [26]. A path MI that is sensitive to temporal effects has been recently introduced for a class of chemical reaction networks [13]. The MI is often interpreted as a measure for the amount of input states that can be resolved accurately [27]. From an engineering perspective maximizing the path MI between sensor and actuator could thus define a design principle to construct cellular circuits. Information theory has been used to reveal fundamental physical limits of information transmission guiding our understanding of gene regulatory motifs as information processing units [12, 28].

The standard transcription model with a two-state promoter. a) A random telegraph promoter switches between its active and its inactive state. Only in the active state, mRNA is synthesized. The promoter stays active, so multiple mRNA strands can be transcribed in one active period. The transcription reaction is of first order with fixed rate. b) Experimentally only the transcription events are observed. The promoter is modelled as a context. An observation model for transcription counts is obtained when we marginalize the joint system over the context. The transcription reaction is of order zero with stochastic rate, i.e., a doubly stochastic Poisson process.

Why is the Poisson channel used as a communication model in cells? Chemical reaction networks are standard Markov approaches in the stochastic modeling of gene regulation in cells. Reactants $X_1, \dots, X_n$ are converted, synthesized and degraded in reactions

$$ R_j: \sum_{i = 1}^n\nu_{ij} X_i \overset{r_j}{\longrightarrow} \sum_{i = 1}^n\eta_{ij} X_i, \quad j = 1, \dots, m. $$

This induces a continuous time Markov chain (CTMC) on the state space $\N_{\geq 0}^n$ of molecule counts $\vec{k} = (k_1, \dots, k_n)$, where the non-diagonal entry $(\vec{k} \to \vec{l})$ in the generator matrix is

$$\sum_{j = 1}^m r_j\prod_{i = 1}^n\binom{k_i}{\nu_{ij}}\mathds{1}(\vec{l}-\vec{k} = \vec{\eta}_j - \vec{\nu}_j).$$

The order of reaction $R_j$ is $\sum_{i = 1}^n \nu_{ij}$. For instance, the standard transcription model with the two-state promoter (Fig. <ref>a) contains a first-order synthesis reaction

$$ P_{\mathrm{on}} \to  P_{\mathrm{on}} + \mathrm{mRNA}.$$

The homogeneous Poisson process corresponds to a zero-order reaction. First-order reactions can be interpreted as zero-order reactions with stochastic modulation, i.e., as Poisson channel with linear $\intensityFct(\msg_\ti)$ (Fig. <ref>b).
More generally, the Poisson channel accounts for the discrete nature of reactions $R_j$ with the reaction counts $\out_j(\ti)$ up to time $\ti$ being the Poisson channel output. Suited for the low copy number regime, it can bring to attention the bottleneck character of sensor molecules that sense a continuous signal but whose synthesis events are restricted to the discrete regime [12, 29]. Adhering to our goal of including temporal effects, the Poisson channel is suited because it is capable of modeling time-varying inputs and outputs.

Computing the path MI over Poissonian channels is difficult. Many works of research resort to basing the MI on single time-point marginals ignoring any encoding in the temporal profile [30, 31, 22].
Other approaches are the following: (i) Gaussian approximations of the input to make use of analytic results [12], (ii) Monte Carlo estimators [32, 33, 13]. (iii) Not rarely, the intractability decoyed researchers to resort to other channels such as the Gaussian channel [34, 35]. We contribute to information theory in biology by presenting a Monte Carlo-free numerical computation of the path mutual information rate of a Poisson channel for Markovian input with low state number.

§.§ Approximate filters and approximate marginal simulation

The study of ACID can be beneficial beyond its appearance in the information rate of the Poisson channel.
Markov-modulated Poisson processes do not only serve as communication channel models. They provide observation models for open systems, i.e., the counting process can be regarded as a subsystem that is embedded in a heat bath [36] or environment [37]. The Markovian environment modulates the intensity of the observed subsystem.

In the biological cell, for instance,  the rate of transcription synthesis events is modulated by multiple factors [38]. Recruitment of polymerase and transcription factors as well as unwinding of the DNA strand, all contribute to the activation of a promoter state, where transcription is initiated [39]. The combination of these factors can be regarded as a random environmental process modulating the rates in the observed main process, here transcription synthesis (Fig. <ref>b).

An observation model of the subsystem requires modeling the unobserved context, e.g., as Markovian environment.
After marginalizing over the environment, the CI describes the main process in an uncoupled way.
Sample trajectories, i.e., Monte Carlo samples, can be obtained by simulating the main process as a self-exciting counting process whose rate is the CI. In practice, this is done via thinning [40] or by the inverse transform method to sample sojourn times.
This marginal simulation replaces the co-simulation of the environment and promises a potential speed-up. Only the effective rate $\CIntensity_\ti$ is used, neglecting fluctuations of the environmental process that are not transmitted to the main process. However, unclosed conditional moment equations pose a serious challenge to achieving the goal of fast marginal simulation [37]. This problem is addressed by approximate filters obtained from conditional moment closure, assumed density filtering, variational inference, entropic matching or projection [41, 42, 3]. Prominently, the Hawkes process can be obtained from projection onto the class of linear estimators and can thus be regarded as an approximate filter. We call the marginal simulation with an approximate filter approximate marginal simulation.

The process that results from approximate marginal simulation is self-exciting, but in general not a Markov-modulated Poisson process any more. Do characteristics of $\out_\ti$, such as the mean or variance evolution, stay invariant under approximate marginal simulation? A full, but often intractable characteristic is the path measure of $\out_\ti$. The ACID provides a partial characteristic that exceeds first- and second-order analysis but is less complex than the path measure. It can assist as a discrimination tool to detect dissimilar approximate filters. We demonstrate this using the Hawkes process.

§ METHOD

Let $\out_\ti$ be a counting process. Its canonical filtration is $\sigmalg_\ti^\out = \sigma(\out_s: s \le \ti)$. Denote by $\CIntensity_\ti$ the $\sigmalg_\ti^\out$-intensity of the process $\out_\ti$. We refer to it as the conditional intensity (CI).

§.§ The backward recurrence time parametrization

Throughout this work we consider a counting process $\out_\ti$ whose conditional intensity $\CIntensity_\ti$ is parametrized by a process $(\intime(\ti), \state(\ti))$ in the following form. The scalar $\intime$ is the time since the latest jump of $\out$, i.e., the backward recurrence time [43].

The second component $\state \in \R^\num$ is a (possibly multi-variate) sufficient statistic, that possesses three properties

(A1) it is constant between jumps, i.e. $\state(\ti) = \state(\ti - \intime(\ti))$

(A2) there is a deterministic $\cintensity$, satisfying $\CIntensity_\ti = \cintensity(\intime(\ti), \state(\ti))$.

(A3) there is a deterministic $\target$, satisfying $\target(\intime(\ti-), \state(\ti-)) = \state(\ti)$ at jump times $\ti$ of $\out$.

We further request $(\intime(\ti), \state(\ti))$ to be ergodic in the following. An example of a non-ergodic process $(\intime(\ti), \state(\ti))$ satisfying the three properties was investigated in [44].
The process $(\intime(\ti), \state(\ti))$, denoted BReT-P, combines three conditions which jointly turn it into a piecewise-deterministic Markov process [45]:

(B1) The evolution equation of $(\intime, \state)$ between jumps reads

\begin{equation}
\dot\intime = 1, \, \dot\state = 0.
\label{eq:process_equation}
\end{equation}

Thus the process evolves deterministicly from the current state when there are no jumps.

(B2) Jumps of $(\intime(\ti), \state(\ti))$ and $\out$ occur simultaneously with intensity $\CIntensity_\ti = \cintensity(\intime(\ti), \state(\ti))$, only depending on the current state, not the history.

(B3) At jumps, the backwards recurrence time is updated to $\intime(\ti) = 0$ and the sufficient statistic targets $ \state(\ti) = \target(\intime(\ti-), \state(\ti-))$. The jump targets only depend on the current state.

Next, we show that the CI of a Markov-modulated Poisson process is of this form. After the joint asymptotic distribution of $(\intime, \state)$ is found (see section <ref>), the ACID is determined via the transform $\cintensity(\intime, \state)$ (see section <ref>).

§.§ Class of Markov-modulated Poisson processes

Let $\msg$ be a continuous time Markov chain (CTMC) with finitely many states in $\msgStateSpace$ and $\intensityFct: \msgStateSpace \to [0, \infty)$ an intensity mapping.
Let $\out$
be a doubly stochastic Poisson process with intensity $\intensity_\ti = \intensity(\msg_\ti)$.

Denote by $\MasterEq$ the generator of $\msg$. Then the filtering distribution $\CProb_\ti(\msgState) := \Prob[\msg_\ti = \msgState\vert \sigmalg^\out_\ti], \msgState \in \msgStateSpace$ evolves as follows, according to Snyder [3],

\begin{equation}
\frac{\dd}{\dd \ti} \CProb_\ti(x) = (\MasterEq \CProb_\ti)(x) -  (\intensityFct(x)-\CIntensity_\ti) \CProb_\ti(x)
\label{eq:between_jumps}
\end{equation}

between jumps and is updated to

\begin{equation}
\CProb_\ti(\msgState) = \frac{\intensityFct(\msgState)\CProb_{\ti-}(\msgState)}{\CIntensity_{\ti-}}
\label{eq:at_jumps}
\end{equation}

if $\out_\ti = \out_{\ti-} + 1$. Taken together, the Snyder filter reads

\begin{equation}
\dd \CProb_\ti(x) = (\MasterEq \CProb_\ti)(x) \dd \ti + \frac{(\intensityFct(x) - \CIntensity_\ti) \CProb_\ti(x) }{\CIntensity_\ti} \{ \dd \out_\ti - \CIntensity_\ti \dd \ti \}.
\label{eq:filter}
\end{equation}

By means of the filtering distribution, the conditional intensity is computed as
$$\CIntensity_\ti= \E [\intensity(\msg_\ti)\vert \sigmalg^\out_\ti] = \sum_\msgState \intensityFct(\msgState)\CProb_{\ti}(\msgState).$$
It is also denoted as the filter mean or causal conditional mean estimate.

We now consider $\state(\ti) := \CProb_{\ti} \in [0,1]^{\vert \msgStateSpace \vert}$ at jump times and $\state(\ti) = \state(\ti - \intime(\ti))$ else, and derive that the process $(\intime(\ti), \state(\ti))$ is a BReT-P of $\CIntensity_\ti$. In order to do so, we introduce auxiliary functions $\CIprojection, \auxODE$ and $\jtarget$ to define the functions $\cintensity$ and $\target$ in this context. Let $\CIprojection(\cprob):=\sum_{\msgState}\intensityFct(\msgState)\cprob(\msgState)$ be the mean functional. Observe that, by the relation $\CIntensity_\ti = \CIprojection(\CProb_{\ti})$, the ODE system (<ref>) is closed and autonomous. Denote by $\intime \mapsto \auxODE(\intime, \cprob)$ its solution  with initialization $\CProb_0 = \cprob \in [0,1]^{\vert \msgStateSpace \vert}$.
The jump update $\jtarget \colon \R^{\vert \msgStateSpace \vert} \to \R^{\vert \msgStateSpace \vert}$ is $\jtarget(\cprob)(x) = \frac{\intensityFct(\msgState)\cprob(\msgState)}{\CIprojection(\cprob)}$. Then $\cintensity := \CIprojection \circ \auxODE$ and $\target := \jtarget \circ \cintensity$ satisfy the properties (A2) and (A3).

It is crucial to note that we can indeed consider $\out_\ti$ to be jumping with intensity $\CIntensity_\ti$ instead of $\intensity_\ti$. In condition (B2), if we replaced $\CIntensity_\ti$ by $\intensityFct(\msg_\ti)$ we would lose the self-contained description of $\CIntensity_\ti$ that even disposes of $\out_\ti$ itself. To this end, in process equation (<ref>), replace $\out_\ti$ by $\tilde{\out}_\ti$, where $\tilde{\out}_\ti$ is a self-exciting counting process with intensity $\CIntensity_\ti$. The processes $\tilde{\out}_\ti$ and $\out_\ti$ are equal in distribution. So instead of $\out_\ti$ we consider $\tilde{\out}_\ti$ in the first place and drop the tilde for convenience.

The fact that the Snyder filter is a piecewise-deterministic Markov process has been observed before [46]. For the particle interpretation of this fact we refer to [18].

§.§ Sufficient state variables of joint Markovian progression

We can abstract the BReT-P in subsection <ref> to a more general process class that falls in our considered framework. Consider state variables $(\SuffStat_1(\ti), \dots, \SuffStat_{n_0}(\ti)) = \SuffStat(\ti) \in \R^{n_0}$ that progress as

\begin{equation}
\dot\SuffStat(\ti) = \SuffDyn(\SuffStat(\ti))
\label{eq:SuffODE}
\end{equation}

with the deterministic dynamics $\SuffDyn$ between jumps of $\SuffStat(\ti)$, and at jumps they are updated to

\begin{equation}
\SuffStat(\ti) = \jtarget(\SuffStat({\ti-}))
\label{eq:SuffUpdate}
\end{equation}

with the update function $\jtarget = (\jtarget_1, \dots, \jtarget_{n_0}) \colon \R^{n_0} \to \R^{n_0}$.
Let further be $\CIprojection \colon \R^{n_0} \to \R$ a deterministic functional of the state variable, such that $\CIntensity_\ti = \CIprojection(\SuffStat(\ti))$. We call $\SuffStat_1, \dots, \SuffStat_{n_0}$ sufficient state variables of joint Markovian progression, because they form a Markov process by arguments analogous to (B1)-(B3) and are sufficient in that the conditional intensity can be computed from them. Denote by $\intime \mapsto \auxODE(\intime, \suffStat^0)$ the solution of eq. (<ref>) with initial value $\SuffStat_0 = \suffStat^0$.
Suppose that $\jtarget_{n+1} \equiv \suffStat^0_{n+1}, \dots, \jtarget_{n_0}\equiv \suffStat^0_{n_0}$ are constant, i.e., the values $\SuffStat_{n+1}, \dots, \SuffStat_{n_0}$ are reset to the same values $\suffStat^0_{n+1}, \dots, \suffStat^0_{n_0}$ at any jump. We now construct the process $(\intime(\ti), \state(\ti))$. As sufficient statistic it suffices to define $\state(\ti) = (\SuffStat_1(\ti), \dots, \SuffStat_n(\ti))$ at jumps and $\state(\ti) = \state(\ti - \intime(\ti))$ else, automatically satisfying (A1). The extension operator $\extension(\state) := (\state, \suffStat^0_{n+1}, \dots, \suffStat^0_{n_0})$ concatenates the current state with the constants. The truncation operator $\updateProj(\suffStat^0) = (\suffStat^0_1, \dots, \suffStat^0_n)$ projects onto the first $n$ components. The functions $\cintensity(\intime, \state) = \CIprojection \circ \auxODE(\intime, \extension(\state)))$ and $\target(\intime, \state) = \updateProj \circ \jtarget(\auxODE(\intime, \extension(\state)))$ satisfy the properties (A2) and (A3).

We return to section <ref> and exclude the state variables with constant reset value. Namely, conservation of probability mass and zero-states reduce the dimension of the sufficient statistic.

* In eq. (<ref>) the evolution equation of the last $\msgState$ can be replaced by the trivial evolution of $\sum_{\msgState} \CProb_\ti(\msgState)$. The value of this sum is constantly $1$. Hence, the number of state variables $\num_0$ is at most $\vert \msgStateSpace \vert -1$. A reparametrization may further decrease the number.

* Call any $\msgState \in \msgStateSpace$ with $\intensityFct(\msgState) = 0$ a zero-state. For zero-states, the reset value of the corresponding $\CProb_\ti(\msgState)$ in eq. (<ref>) is $0$, i.e., $\jtarget$ is constant for these components. Hence, the conditional probabilities for the zero-states need not be tracked in the sufficient statistic $\state$. The dimension of $\state$ can be reduced to

\begin{equation}
\num \le \vert \{\msgStateSpace  \colon \intensityFct(\msgState) > 0\} \vert -1.
\label{eq:zero-states}
\end{equation}

§.§ BReT-P examples

Model    Type    $\SuffStat(\ti)$     $\num_0$     $\num$     $\CIprojection$     $\auxODE$     $\SuffDyn$     $\jtarget$

Random telegraph   MM    $\CProb_\ti(1)$    $1$   $0$    $\gain\cprob_1$         Eq. (<ref>)     $1$

Dark current   MM    $\CIntensity_\ti$    $1$   $1$    $ \id$       Eq. (<ref>)   $\CIntensity^{-1}(\CIntensity - \leakage)(\ONState - \CIntensity) + \CIntensity$

Double On   MM    $(\DonsoffFrac_\ti,\condmean_\ti, \CProb_\ti(\msgState_3),$
$\dots, \CProb_\ti(\msgState_{\vert \msgStateSpace \vert -1}))$      $\vert \msgStateSpace \vert - 1$    $1$    $ \gain \z$    -   Ito from
Eq. (<ref>)   $(\Donsofffrac, 1, 0, \dots, 0)$

Double On Single Off   MM   $(\DonsoffFrac_\ti,\condmean_\ti)$    $2$   $1$    $\gain \z$    -   Eq.(<ref>)-(<ref>)   $(\Donsofffrac,1)$

Hawkes   SE    $\CIntensity_\ti$    $1$   $1$    $\id$       Eq. (<ref>)   $\CIntensity + \HawkesJump$

Gamma filter   SE    $(\GammaMean_\ti, \GammaVar_\ti)$    $2$   $2$    $\gain \gammamean$    -    Eq. (<ref>)   $(\gammamean, \Gammavar) + (\Gammavar\gammamean^{-1}, \Gammavar^2\gammamean^{-2})$

The type indicates whether the model is a Markov-modulated (MM) Poisson process or self-exciting (SE) counting process. The minus ($-$) in the column $\auxODE$ indicates that no analytic solution is available for $\intime \mapsto \auxODE(\intime, \suffStat^0)$.

We list examples of parametrizations in order to convey the intuition behind the sufficient state variables of joint Markovian progression. The examples introduced here serve as our case studies in section <ref>, highlighting different aspects of the BReT-P method and usage of ACID. In case studies 1) to 3) ACID is used to evaluate the mutual information rate via (<ref>). Examples 4) and 5) show how the analysis of ACID exceeds first and second order analysis, both on the CI level and the level of $\out_\ti$. Contrasting examples 2) and 4) illustrates how ACID can discriminate approximate filters. The table <ref> summarizes the computational details of the examples.

§.§.§ Random telegraph model

The input $\msg$ is a random telegraph model on $\msgStateSpace := \{0,1 \}$ with On and Off rates $\LeakageOn, \LeakageOff$ and $\intensityFct(\msgState) = \gain\msgState$. By remark <ref>(i) $\CProb_\ti(1)$ is a sufficient state variable (i.e., $\num_0 = 1$). It evolves according to eq. (<ref>) as

\begin{equation}
\dot \CProb_\ti(1) = \LeakageOn - (\LeakageOff + \LeakageOn + \gain) \CProb_\ti(1) + \gain\CProb_\ti(1)^2.
\label{eq:logistic}
\end{equation}

For the conditional intensity, we get

$$ \CIntensity_\ti = c \cdot \CProb_\ti(1) + 0\cdot (1-\CProb_\ti(1))$$

which means $\CIprojection(\cprob_1) = c\cprob_1$. The zero-state reduces the dimension. Indeed, by eq. (<ref>) the reset is $\jtarget(\cprob_1) = 1$, so the dimension of $\state$ can be chosen $\num = 0$ as discussed in subsection <ref>, compare also eq. (<ref>). Consequently, the scalar $\intime(\ti)$ suffices as BReT-P. Let $\intime \mapsto \auxODE(\intime)$ be the solution of the ODE with initial value $\cprob_1 = 1$, then $\CIntensity_\ti = \cintensity(\intime(\ti)) = \CIprojection \circ \auxODE(\intime(\ti)) = \gain \auxODE(\intime(\ti))$.

The BReT-P method helps in analyzing how $\IR(\msg, \out)$ depends on the system parameters $\LeakageOn, \LeakageOff, \gain$, see paragraph <ref>, with the goal to answer questions of optimality under constraints $0 < \LeakageOn \le \OnConstraint, 0 < \LeakageOff \le \OffConstraint$.

§.§.§ Random telegraph model with dark current

Let $\msg$ be a random telegraph model as before, but this time, $\intensityFct(0) =:\leakage > 0 $ introduces a dark current. Define $\ONState := \intensityFct(1) > \intensityFct(0)$ and the amplitude $\dynamicRange := \ONState - \leakage$. By remark <ref>(i) again $\num_0 = 1$. Choose $\CIntensity_\ti = \leakage +  \CProb_\ti(1) \dynamicRange$ as sufficient state variable of Markovian progression. Being an affine-linear transform of eq. (<ref>), it evolves as

\begin{equation} \begin{split}
\dd \CIntensity_\ti = \left\{\LeakageOn\dynamicRange - (\LeakageOff + \LeakageOn + \dynamicRange)(\CIntensity_\ti - \leakage) \right.\\
\left.+ (\CIntensity_\ti - \leakage)^2\right\} \dd \ti + \frac{(\CIntensity_{\ti-} - \leakage)(\ONState-\CIntensity_{\ti-})}{\CIntensity_{\ti-}} \dd \out_\ti.
\end{split}
\label{eq:Leakage}
\end{equation}

(By their affine-linear relation, $\CIntensity_\ti$ and $\CProb_\ti(1)$ are equivalent choices.)
From eq. (<ref>) we can only conclude $\num \le 1$. Thus, we have $(\intime(\ti), \state(\ti))$ with $\state(\ti) = \CIntensity_\ti$ at jumps as BReT-P.
BReT-P helps in investigating how dark current alters $\IR(\msg, \out)$, see paragraph <ref>.

§.§.§ CTMC with two On states
Let $\msg$ be an arbitrary ergodic CTMC with $n$ states and $\intensityFct: \msgStateSpace \to \{ 0,\gain \}$ be binary with two active, i.e., non-zero states $\msgState_1, \msgState_2$. We refer to the model as double On (DOn). Then introduce the conditional probability of being in an active state

$$\condmean_\ti = \CProb_\ti(\msgState_1) + \CProb_\ti(\msgState_2)$$
and $\DonsoffFrac_\ti:= \CProb_\ti(\msgState_1)/\condmean_\ti$ the contribution of $\msgState_1$ to this conditional probability. Then $\SuffStat_\ti := (\DonsoffFrac_\ti, \condmean_\ti, \CProb_\ti(\msgState_3), \dots, \CProb_\ti(\msgState_{\vert \msgStateSpace \vert - 1}))$ are sufficient variables of joint Markovian progression. The reset value at jumps is $(\DonsoffFrac_{\ti-}, 1, 0, \dots, 0)$. Since all but the first component are set to a constant at jumps, the sufficient statistic $\state(\ti) = \DonsoffFrac_\ti$ at jumps has dimension $\num = 1$. The progression of  $\SuffStat_\ti$ can be found by the chain rule from eq. (<ref>). We only elaborate on it for the circular CTMC with double On state and single Off state (DOnSOff), depicted in figure <ref>.

[->, >=latex', auto, semithick, node distance=1/3*]
every state=[fill=white,draw=black,thick,text=black,scale=1]
[state]    (A)                     Off;
[state]    (B)[right of=A]   On;
[state]    (C)[right of=B]   On;
(A) edge     node$\DonsoffOn$     (B)
(B)    edge    node$\DonsoffR$      (C)
(C)   edge[bend left, below]     node$\DonsoffOff$     (A);

State diagram for the Double On Single Off (DOnSOff) model. The three-state model is Markovian. A refractory second active state realizes a non-exponential sojourn time in On. For $\DonsoffR = \DonsoffOff$ the sojourn time is an Erlang distribution. The Off sojourn time remains exponential.

Let $\DonsoffOn, \DonsoffR, \DonsoffOff$ be the transition rates of going from inactive to first active, first to second active and second active back to inactive, then the generator matrix is of the form

\begin{equation}
\MasterEq = \begin{pmatrix}
-\DonsoffR & \DonsoffR& 0\\
0 & -\DonsoffOff & \DonsoffOff\\
\DonsoffOn & 0 & -\DonsoffOn
\end{pmatrix}.
\end{equation}

The evolution of the conditional probabilities $\CProb_\ti(\msgState_1)$ and $\CProb_\ti(\msgState_2)$ is

\begin{align}
&\dot\CProb_\ti(\msgState_1) \nonumber\\ &= \DonsoffOn(1 - \condmean_\ti) - \DonsoffR\CProb_\ti(\msgState_1) - \gain(1 - \condmean_\ti))\CProb_\ti(\msgState_1)\\
&\dot\CProb_\ti(\msgState_2)\nonumber\\ &= \DonsoffR\CProb_\ti(\msgState_1) - \DonsoffOff\CProb_\ti(\msgState_2) - \gain(1 - \condmean_\ti))\CProb_\ti(\msgState_2).
\end{align}

Via the function $\varphi(\cprob_1, \cprob_2) = (\cprob_1/(\cprob_1 + \cprob_2), \cprob_1 + \cprob_2) = (\Donsofffrac, \z)$ this transforms to

\begin{align}
\dot \DonsoffFrac &= - \DonsoffR \DonsoffFrac + \DonsoffOn \frac{(1-\condmean)(1-\DonsoffFrac)}{\condmean} + \DonsoffOff \DonsoffFrac(1-\DonsoffFrac) \label{eq:Donsoff_frac}\\
\dot \condmean&= - \gain (1-\condmean) \condmean + \DonsoffOn (1-\condmean) - \DonsoffOff (1-\DonsoffFrac) \condmean. \label{eq:Donsoff_cintensity}\
\end{align}

The update function according to eq. (<ref>) is $\jtarget(\Donsofffrac, \z) = (\Donsofffrac, 1)$, i.e., constant in the second component. Figure <ref> visualizes how trajectories for different initial values $[\DonsoffFrac(\intime), \condmean(\intime)]_{\intime = 0} = [\state, 1]$ evolve. The coordinate system shows the plane $(\cprob_1, \cprob_2)$, while the grid indicates the transformed radial-like coordinates $(\Donsofffrac, \z)$.

The Double On Single Off input. Trajectories $(\CProb_\intime(\msgState_1), \CProb_\intime(\msgState_2))$ evolve in the $(\cprob_1, \cprob_2)$-plane with $\cprob_1+ \cprob_2 \le 1$. The grey grid indicates the change of coordinates $(\Donsofffrac, \z) = (\cprob_1/(\cprob_1+ \cprob_2), \cprob_1+ \cprob_2)\in [0,1]^2$. The initial value is $[\Donsofffrac_0, \z_0] = [\state, 1]$. Dots indicate the temporal evolution, i.e. trajectory values at equally spaced time points. The equilibrium is marked by the cross. Parameter values were $\DonsoffOn = 0.4, \DonsoffR =\DonsoffOff = 0.8, \gain = 1, \state = 0.2, 0.4, 0.6$

The DOnSOff model is equivalent to a binary semi-Markov process with exponential sojourn time in the Off state, while the sojourn time in the On state is the convolution of two exponential distributions. It serves as an example of a non-Markovian binary input to the Poisson channel whose $\IR(\msg, \out)$ is compared to the Markov case, see case study $\ref{cs:donsoff}$.

§.§.§ Hawkes process
The Hawkes process [19] is a self-exciting counting process whose conditional intensity evolves like

\begin{equation}
\hat\intensity_\ti = \HawkesBase + \int_0^\ti \HawkesKernel(s, \ti) \dd \out_s,
\label{eq:general_Hawkes}
\end{equation}

where $\HawkesBase \geq 0$ is a base intensity. For the standard Hawkes process with exponential kernel $\HawkesKernel(s,\ti) = \HawkesJump e^{-\HawkesDecay(\ti - s)}$, the CI turns out to be of Markovian progression [20, 17]

\begin{equation}
\dd \hat\intensity_\ti = -\HawkesDecay(\hat\intensity_\ti - \HawkesBase) \dd \ti + \HawkesJump \dd \out_\ti.
\label{eq:Hawkes}
\end{equation}

The sufficient statistic $\state(\ti)$ can be chosen as $\CIntensity_\ti$ at jump times. The functions $\cintensity(\intime, \state) = \HawkesBase + e^{-\HawkesDecay \intime} (\state - \HawkesBase)$ and $\target(\intime, \state) = \cintensity(\intime, \state) + \HawkesJump$ satisfy (A2) and (A3). The equilibrium distribution of $\CIntensity_\ti$, i.e., the ACID, is mentioned in the literature [20, 17], partly because it helps in realizing a stationary Hawkes process. A Hawkes process whose initial intensity value is drawn from ACID is in the stationary regime from the start. We demonstrate how ACID contains more information than the mean and variance of $\CIntensity_\infty$, see case study <ref>. Oakes [20] suggested that the ACID can be found by iteratively solving an integral equation, but does not elaborate more on it. Daley & Vere-Jones <cit.> provide a manual for setting up the integral equation and solving it with the method of steps. We are not sure whether Oakes meant this method by "iterative solution". We provide a different approach using BReT-P and a fixed point iteration.

Besides being an interesting counting process model on its own that has seen a variety of applications in different fields of research, the Hawkes process can also be regarded as an approximate filter. It is obtained by projection onto the class of linear estimators under the quadratic criterion and we therefore also call it the optimal linear estimator. In case study <ref> we elaborate on this and employ the ACID to discriminate between the exact Snyder filter and the optimal linear filter in <ref>.

§.§.§ The Gamma filter
The Gamma filter [37, 41] is an approximate filter obtained from conditional moment closure. It departs from the linear modulation by a Cox-Ingersoll-Ross process [47] with stationary CIR-mean $\GammaMu$ and CIR-autocovariance function $\GammaSigma e^{-\GammaDecay \ti}$, i.e., $\msg_\ti \sim \mathrm{CIR}(\GammaMu, \GammaSigma, \GammaDecay)$ and $\out_\ti$ is a doubly stochastic Poisson process with $\intensity_\ti = \gain \msg_\ti$. The Gamma filter is the assumed density filter, where $\msg_\ti\vert \out_{[0,\ti]}$ is assumed to be Gamma distributed. With two degrees of freedom for the Gamma distribution, two equations sufficiently describe the filter. One governs the mean, one governs the variance.
Expressing the Gamma's third centered moment in terms of mean and variance justifies the replacement $\E[(\msg_\ti - \GammaMu)^3 \vert \sigmalg^\out_\ti] = \frac{2\E[(\msg_\ti - \GammaMu)^2 \vert \sigmalg^\out_\ti]^2}{\E[\msg_\ti \vert \sigmalg^\out_\ti]}$. Consider an approximate marginal simulation of $\out_\ti$ that uses the Gamma filter. This yields a self-exciting process $\out_\ti$ with $(\GammaMean_\ti, \GammaVar_\ti)$ mimicking $(\E[\msg_\ti \vert \sigmalg^\out_\ti], \E[(\msg_\ti - \GammaMu)^2 \vert \sigmalg^\out_\ti])$, the CI $\CIntensity_\ti = \gain \GammaMean_\ti$ and

\begin{align}
\dd\GammaMean_\ti &= \{-\GammaDecay (\GammaMean_\ti - \GammaMu) - \gain \GammaVar_\ti \} \dd \ti + \frac{\GammaVar_{\ti-}}{\GammaMean_{\ti-}}\dd \out_\ti \nonumber\\
\dd\GammaVar_\ti &= \{-2\GammaDecay (\GammaVar_\ti - \frac{\GammaSigma}{\GammaMu}\GammaMean_\ti) - 2\gain\frac{\GammaVar_\ti^2}{\GammaMean_\ti} \} dt + \frac{\GammaVar_{\ti-}^2}{\GammaMean_{\ti-}^2}\dd \out_\ti.
\label{eq:Gamma_filter}
\end{align}

It is instructive to contrast approximate marginal simulation of $\out_\ti$ with estimation of $\msg_\ti$. Both perspectives can make use of the equations (<ref>). In estimation the process $\out_\ti$ is self-exciting with intractable CI $\gain \E[\msg_\ti\vert \sigmalg^\out_\ti]$, while in approximate marginal simulation $\out_\ti$ is by definition self-exciting with $\gain \GammaMean_\ti$.
We proceed with the latter perspective.
While the CI $\CIntensity_\ti$ alone is not of Markovian progression, the joint $\SuffStat_\ti = (\GammaMean_\ti, \GammaVar_\ti)$ is, yielding $\num = \num_0 = 2$.

With $\state(\ti) = (\GammaMean_\ti, \GammaVar_\ti)$ at jump times, the sufficient statistic can be defined.

The Gamma filter and optimal linear filter are obtained from different closure schemes. In case study <ref> the ACID is employed to discriminate between approximate marginal simulations with either one. The case study of the Gamma filter indicates the limitations of ACID as a discrimination tool for approximate filters. It, however, informs a decision, in which parameter regime to replace the Gamma by the optimal linear filter when approximately computing the information rate.

The examples illustrate how the concept of CI unifies the two pictures:

* the conditional mean estimate in a Markov-modulated Poisson process and

* the history-dependent intensity in a self-exciting process, e.g. Hawkes process and marginal simulation via approximate filters.

§.§ The master equation and stationarity condition

We now aim for the asymptotic distribution of $(\intime(\ti), \state(\ti))$.
The master equation for the piecewise-deterministic Markov process $(\intime(\ti), \state(\ti))$, derived from the process equation $\eqref{eq:process_equation}$, reads [21]

\begin{equation}
\partial_\ti \density(\ti, \intime, \state) = - \partial_{\intime} \density(\ti, \intime, \state) -  \cintensity(\intime, \state) \density(\ti, \intime, \state).
\label{eq:Master_equation}
\end{equation}

At stationarity the left-hand side vanishes, yielding

\begin{equation}
\partial_{\intime} \density(\intime, \state) = - \cintensity( \intime, \state) \density(\intime, \state).
\label{eq:stationarity}
\end{equation}

Let $\Density$ solve this equation for $\Density(0, \state) = 1$. Since eq. (<ref>) is linear, the true density then satisfies

\begin{equation}
\density(\intime, \state) = \density(0, \state)\Density(\intime, \state).
\label{eq:factorization}
\end{equation}

§.§ Integral boundary condition and normalization condition

The true $\density(0,\state) =: \idensity(\state)$ has to satisfy an influx boundary condition. We state two versions of the condition in the form of integral equations.

For any $\stateSet \in \mathcal B(\R^n)$ it holds

\begin{equation}
\int_\stateSet \idensity(\state) \dd \state = \int_{\target(\intime, \state) \in \stateSet} \cintensity(\intime, \state)\Density(\intime, \state) \idensity(\state) \dd \state \dd \intime.
\label{eq:version1}
\end{equation}

Let $t >h> 0$. The probability $\Prob[\intime(\ti)\in [0,h), \state(\ti) \in \stateSet]$ can be written in two ways (up to order $o(h)$), first

$$ \int_0^h \int_\stateSet \density(\intime, \state, \ti) \dd \state \dd \intime $$

and second, since we know a jump must have occurred in $(\ti-h, \ti]$

\begin{align*}
&\int_{\ti - h}^\ti \int \density(\text{jump at }h'\vert \state(h') = \state', \intime(h') = \intime') \times\\
&\qquad \qquad \density(\intime, \state, h') \Id(\target(\intime', \state') \in \stateSet) \dd \state' \dd \intime'\dd h' \\
&= \int_{\ti - h}^\ti \int\cintensity(\intime', \state') \density(\intime', \state', h') \Id(\target(\intime', \state') \in \stateSet) \dd \state' \dd \intime'\dd h'.
\end{align*}

Then dividing by $h$ and letting $h \to 0$ gives the equality

\begin{equation*}
\int_\stateSet \density(0,\state, \ti) \dd \state = \int_{\target(\intime, \state) \in \stateSet} \cintensity(\intime, \state)\density(\intime,\state, \ti) \dd \state \dd \intime.
\end{equation*}

If we drop the $\ti$ because of stationarity and use eq. (<ref>) on the right-hand side, we get the result.

The differential version of this for the special case $n = 1$ is

For each $\stateto, \statefrom$ let $\intime_i(\stateto, \statefrom)$ for $i = 1, \dots, N(\stateto, \statefrom)$ be an enumeration of the solutions to $\target(\intime, \statefrom) = \stateto$. Assume $\partial_\intime \target(\intime_i(\stateto, \statefrom), \statefrom) \neq 0$ for all $i$.
Here, $N(\stateto, \statefrom)$ is the number of such solutions. Then it holds

\begin{equation}
\idensity(\stateto) = \left.\int \sum_{i = 1}^{N(\stateto,\statefrom)} \frac{\cintensity(\intime, \statefrom)\Density(\intime, \statefrom)}{\vert \frac{\partial}{\partial \intime} \target(\intime, \statefrom) \vert} \idensity(\statefrom) \right\vert_{\intime = \intime_i(\stateto, \statefrom)} \dd \statefrom .
\label{eq:version2}
\end{equation}

Split $\{1, \dots,N(\stateto, \statefrom)\}$ into $\posSet$ and $\negSet$ depending on whether $\intime \mapsto \target(\intime, \statefrom) - \stateto$ has a sign change from $-$ to $+$ or from $+$ to $-$ at $\intime_i(\stateto, \statefrom)$. We choose $\stateSet = (-\infty, \state]$ in eq. (<ref>) and take the derivative $\partial_\state$ on both sides. Then on the left-hand side we get $\idensity(\state)$. For the right-hand side define $\integrand(\intime, \state) := \cintensity(\intime, \state)\Density(\intime, \state) \idensity(\state)$. There exists some choice $\posSet = \{j_1, \dots, j_{\num_+(\stateto, \statefrom)} \}$ or $\posSet = \{j_1, \dots, j_{\num_+(\stateto, \statefrom)-1} \}, \infty = j_{\num_+(\stateto, \statefrom)}$ and $\negSet = \{i_1, \dots, i_{\num_+(\stateto, \statefrom)} \}$ or $\negSet = \{i_2, \dots, i_{\num_+(\stateto, \statefrom)}\}, i_1 = 0$ for which the right-hand side can be computed as

\begin{align*}
&\partial_\state\int \sum_{k = 1}^{\num_+(\stateto, \statefrom)}\int_{\intime_{i_k}(\stateto, \statefrom)}^{\intime_{j_k}(\stateto, \statefrom)} \integrand(\intime, \statefrom) \dd \intime \dd \statefrom\\
&=\int \sum_{k = 1}^{\num_+(\stateto, \statefrom)} - \partial_\state \intime_{i_k}(\stateto, \statefrom) \integrand(\intime_{i_k}(\stateto, \statefrom), \statefrom) \\&\qquad \qquad \quad+ \partial_\state \intime_{j_k}(\stateto, \statefrom) \integrand(\intime_{j_k}(\stateto, \statefrom), \statefrom) \dd \statefrom\\
&=\int \sum_{i \in \posSet \cup \negSet} \sgn\{\partial_\intime \target(\intime_i(\stateto, \statefrom), \statefrom)\} \partial_\state \intime_{i}(\stateto, \statefrom) \\&\qquad \qquad \qquad \qquad \times\integrand(\intime_{i}(\stateto, \statefrom), \statefrom)\dd \statefrom,
\end{align*}

where we used, that the $\partial_\state$-derivatives vanish for the lower and upper limits $\intime_{i_1}(\stateto, \statefrom) = 0$ and $\intime_{j_{\num_+(\stateto, \statefrom)}}(\stateto, \statefrom) = \infty$.

\begin{align*}
1 &= \partial_\stateto\target(\intime_i(\stateto, \statefrom), \statefrom) \\&= \partial_\intime \target(\intime_i(\stateto, \statefrom), \statefrom) \cdot \partial_\stateto\intime_i(\stateto, \statefrom)\\
&= \vert \partial_\intime \target(\intime_i(\stateto, \statefrom), \statefrom)\vert \cdot \sgn\{\partial_\intime \target(\intime_i(\stateto, \statefrom), \statefrom)\}\partial_\stateto\intime_i(\stateto, \statefrom)
\end{align*}

the result follows.

Note that the integral equation might become singular [48] if there exist $\intime_0, \stateto_0, \statefrom_0$, such that $\target (\intime_0, \statefrom_0) = \stateto_0$, but $\partial_\intime \target (\intime_0,\statefrom_0) = 0$. This may occur both at $\intime_0$ with sign changes of $\intime \mapsto \target (\intime_0, \statefrom_0) - \stateto_0$ and without sign changes. In the first case, the equation (<ref>) ought to be interpreted as having an arbitrary value under the integral for this $(\stateto_0, \statefrom_0)$-pair. In the second case, the value under the integral in eq. (<ref>) is defined for every $(\stateto, \statefrom)$. However, for $\statefrom$ close to $\statefrom_0$ the derivative $\partial_\intime \target(\intime_i(\stateto_0, \statefrom), \statefrom)$ will approach $0$ and hence a singularity will appear nonetheless at $(\stateto_0, \statefrom_0)$. The existence of a solution for this singular integral equation needs to be carefully checked.

Theorem <ref> can be generalized to more than one dimension.
Let $\path_i = (\path_{i,1}, \path_{i,2}) \colon (\low_i,\up_i) \times \stateSupport_i \to [0,\infty) \times \stateSupport, i = 1, \dots, \numPath$ satisfy $\target(\path_i(t, \stateto)) = \stateto$ with each $\path_i$ injective and differentiable and $[0,\infty) \times \stateSupport = \sqcup_{i = 1}^\numPath \path_i((\low_i,\up_i) \times \stateSupport_i)$ up to a set of Lebesgue measure zero. Then it holds that

\begin{align}
\idensity(\stateto)  = \sum_{i\colon \stateto \in \stateSupport_i } \int_{\low_i}^{\up_i} &\vert \det D\path_i(t, \stateto) \vert \cintensity(\path_i(t, \stateto)) \times \nonumber\\
&\Density(\path_i(t, \stateto))\idensity(\path_{i,2}(t, \stateto)) \dd t.
\end{align}

The normalization of $\idensity$ is dictated by the mean intensity, i.e.,

\begin{equation}
\int \idensity(\state) \dd \state = \E[\CIntensity_\infty].
\label{eq:normalization_constant}
\end{equation}

Using $\lim_{\intime \to \infty}\density(\intime, \state) = 0$ and eq. (<ref>), we compute

\begin{align*}
\int \idensity(\state) \dd  \state &= \int \density(0, \state) - \lim_{\intime \to \infty}\density(\intime, \state) \dd  \state\\
& = \int \int_0^\infty -\partial_\intime \density(\intime, \state) \dd \intime \dd \state\\
&= \int\int_0^\infty \cintensity(\intime, \state) \density(\intime, \state) \dd \intime \dd \state
\end{align*}

and recognize the right hand side as $\E[\CIntensity_\infty]$.

We emphasize that the strong advantage of the BReT-P lies in aligning all probability inflow terms in the master equation at $\intime = 0$. Our method provides the asymptotic distribution of $(\intime, \theta)$. If the sufficient state variable $\SuffStat$ itself is chosen as the domain variable of the asymptotic distribution, then the stationarity equation (<ref>) grows wider: influx terms will be needed for any values of $\SuffStat$ that jumps can anticipate. The equation will assume a difference-differential form, e.g. <cit.> or <cit.>. The BReT-P circumvents this difference-differential formulation for which solution techniques, like method of steps, would be needed in general. A direct solution technique for one-dimensional $\SuffStat$ can handle the difference-differential formulation. It uses a fixed point treatment, similar to eq. (<ref>). We refer the reader to paragraph <ref> and continue with the BReT-P method here. The appealing simplicity of eq. (<ref>) - being autonomous linear ODEs - comes at the cost of an integral boundary condition ($\intime = 0$) which is more involved.

A further advantage of the BReT-P lies in dimension reduction: The pair $(\intime, \state)$ uniquely informs all state variables $\SuffStat_{\num +1}, \dots, \SuffStat_{\num_0}$ with constant reset values, i.e. via $(\intime, \state) \mapsto \auxODE(\intime, \extension(\state))$. As discussed, we may consequently dismiss $\SuffStat_{\num +1}, \dots, \SuffStat_{\num_0}$ in the sufficient statistic $\state$. Eq. (<ref>) - (<ref>) show that the state variables of Markovian progression $\SuffStat_{1}, \dots, \SuffStat_{\num_0}$ also form a piecewise-deterministic Markov process. In case we apply the Markov theory onto $\SuffStat_{1}, \dots, \SuffStat_{\num_0}$ instead of $(\intime, \state)$, can we dismiss the state variables with constant reset value as well? Generally not, because the mapping $(\SuffStat_{1}(\ti), \dots,\SuffStat_{\num_0}(\ti)) \mapsto (\SuffStat_{1}(\ti), \dots,\SuffStat_{\num}(\ti))$ can be non-injective. In that case there exists no unique mapping $(\SuffStat_{1}(\ti), \dots,\SuffStat_{\num}(\ti)) \mapsto(\SuffStat_{1}(\ti), \dots,\SuffStat_{\num_0}(\ti))$. As an example, see fig. <ref> with $\SuffStat_1 = \DonsoffFrac, \SuffStat_2 = \condmean$. The trajectory starting at $\DonsoffFrac(0) = 0.6$ intersects $\Donsofffrac = 0.5$ twice.

§.§ Discretization

We discretize the integral boundary condition (<ref>). To this end, we assume that $\idensity(\state)$ is supported on $\stateSupport \subseteq \R^\num$. We choose a partition $(\statePartition_\ipart)_{\ipart = 1, \dots, \inum}$ with equivalent volume $\vol(\statePartition_\ipart) = \vol(\statePartition_\iPart) =: \weight$. For $\num = 1$, this is an equidistant partition. We discretize $\idensity(\state)$ as

\begin{equation}
\idensity(\state) = \sum_{\ipart =1}^\inum \idensityCoeff_\ipart \Id_{\statePartition_\ipart}(\state)
\end{equation}

with unknowns $\idensityCoeff_\ipart$ and discretize $\cintensity(\intime, \state)$, respectively $\Density(\intime, \state)$, as

\begin{equation}
\cintensity(\intime, \state) = \sum_{\ipart =1}^\inum \cintensity(\intime, \state_\ipart) \Id_{\statePartition_\ipart}(\state) , \quad \Density(\intime, \state) = \sum_{\ipart =1}^\inum \Density(\intime, \state_\ipart) \Id_{\statePartition_\ipart}(\state)
\label{eq:discretization}
\end{equation}

for a choice of representatives $\state_\ipart \in \statePartition_\ipart$, e.g. the center of $\statePartition_\ipart$.
Define the border crossing time points $\intime_0(\statefrom) := 0$ and recursively
$\intime_\itime(\statefrom) := \min\{ \intime > \intime_{\itime - 1} : \target(\intime , \statefrom) \in \partial \statePartition_\ipart \text{ for some } \ipart \}$.
Define $\itrap(\itime, \statefrom) :=  \ipart$ if $\target(\intime, \statefrom) \in \statePartition_\ipart$ for $\intime_{\itime-1}(\statefrom) < \intime < \intime_\itime(\statefrom)$. Then for $\stateSet = \statePartition_\ipart$ the equation (<ref>) reads

\begin{align}
\idensityCoeff_\ipart &= \sum_{\iPart = 1}^\inum \left(\sum_{\itime: \itrap(\itime,\statefrom_\iPart) =\ipart }\int_{\intime_{\itime-1}(\statefrom_\iPart)}^{\intime_{\itime}(\statefrom_\iPart)} \cintensity(\intime, \statefrom_\iPart) \Density(\intime, \statefrom_\iPart) \dd \intime \right) \idensityCoeff_\iPart \nonumber\\
&= \sum_{\iPart = 1}^\inum \left(\sum_{\itime: \itrap(\itime,\statefrom_\iPart) =\ipart } \Density(\intime_{\itime-1}(\statefrom_\iPart), \statefrom_\iPart ) - \Density(\intime_{\itime}(\statefrom_\iPart), \statefrom_\iPart ) \right) \idensityCoeff_\iPart.
\label{eq:matrix_version1}
\end{align}

Defining the bracket term as $\OpI_{\ipart, \iPart}$, the equation can be written $\vec\idensityCoeff = \OpI \vec\idensityCoeff$ in matrix form. Observe that by a telescope sum argument and $\Density(0,\statefrom_\iPart) = 1$ the matrix $\OpI$ is left stochastic. If $\OpI$ is quasi-positive, the fixed point equation has a unique non-zero solution by the Perror-Frobenius theorem. Then $\vec\idensityCoeff$ can be approximated by taking a column of a large enough power $\OpI^{2^L}$ and multiplying by the product of weight $\weight$ and normalization constant $\E[\CIntensity_\infty]$.

§.§ The mutual information rate

We assume that $\out_\ti$ is Markov-modulated with intensity $\intensity_\ti$ and CI $\CIntensity_\ti$ given by a BReT-P, e.g., as in section <ref> or <ref>. Let $\intensity_\ti \to \intensity_\infty$ in distribution. The mutual information rate is

\begin{equation}
\IR(\msg, \out) = \E[\phi(\intensity_\infty)] - \int \int_0^\infty \phi(\cintensity(\intime, \state)) \density(\intime, \state) \dd \intime \dd \state.
\label{eq:IR}
\end{equation}

The first term is a finite sum. The outer integral can be approximated by

\begin{equation}
\sum_{\ipart = 1}^\inum \idensity(\state_\ipart) \weight \int_0^\infty  \phi(\cintensity(\intime, \state_\ipart)) \Density(\intime, \state_\ipart) \dd \intime.
\label{eq:IR_discretization}
\end{equation}

Define the partial integral $\cummIR(\ti, \statefrom)$ by

\begin{equation}
\cummIR(\ti, \statefrom) := \int_0^\ti \phi(\cintensity(\intime, \statefrom)) \Density(\intime, \statefrom) \dd \intime.
\end{equation}

In order to proceed from there, we assume that the BReT-P for $\CIntensity_\ti$ is as described in subsection <ref>.
Then for each $\ipart$, the integral in eq. (<ref>) is efficiently solved by the joint ODE system, the dot denoting the $\intime$-derivative,

\begin{align}
\dot \Density(\intime, \state_\ipart) &= - \CIprojection(\suffStat(\intime)) \Density(\intime, \state_\ipart)\label{eq:Density}\\
\dot \suffStat(\intime) &= \SuffDyn(\suffStat(\intime))\\
\dot \cummIR(\intime, \state_\ipart) &= \phi(\CIprojection(\suffStat(\intime)))\Density(\intime, \state_\ipart).\label{eq:cummIR}
\end{align}

The initial conditions are $\Density(0,\state_\ipart) = 1$, $ \suffStat(0) = (\state_\ipart, \suffStat_{\num + 1}^0, \dots, \suffStat_{\num_0}^0)$ and $\cummIR(0,\state_\ipart) = 0$.

If $\cintensity(\intime, \state) \geq \cintensity_0 > 0$ uniformly over $\intime, \state$, then $\Density(\intime, \state) \le e^{-\cintensity_0\intime}$. Let furthermore $\cintensity(\intime, \state) \le R(\state)$. Then

\begin{align*}
\vert \cummIR(\infty, \statefrom) - \cummIR(\ti, \statefrom) \vert &\le \int_\ti^\infty \vert \phi(\cintensity(\intime, \statefrom)) \vert e^{-\cintensity_0 \intime} \dd \intime\\
&= \frac{\max\{\phi(R(\statefrom)), \frac{1}{e}\}}{\cintensity_0} e^{-\cintensity_0 \ti}.
\end{align*}

In order to obtain exponential convergence of $\int \idensity(\statefrom) \cummIR(\ti, \statefrom) \dd \statefrom$ to the integral in eq. (<ref>) for $\ti \to \infty$, it needs to holds that

$$ \int \idensity(\statefrom)\max\{\phi(R(\statefrom)) , \frac{1}{e}\} \dd \statefrom < \infty.$$

Sufficient conditions are, for instance, if $R$ is bounded or if $\state \mapsto R(\state)$ is linear in $\state$ and the second moment of $\idensity$ exists.

As a summary, a numerical approximation of the mutual information rate can be obtained by solving for $\cummIR(\infty, \state_\ipart)$ via eq. (<ref>) - (<ref>) and $\idensity(\state_\ipart)$ via eq. (<ref>). The weights in eq. (<ref>) can be determined on the fly while solving eq. (<ref>) - (<ref>).

§.§ Numerical approximation of the ACID

The ACID is obtained as the distribution $\acid$ of $\cintensity(\intime, \state)$, where $(\intime, \state)$ is distributed according to the density $\density(\intime, \state)$.
We discretize $(0,\infty)$ with mesh size $\Delta \cintensity$ and compute the weights

\begin{align*}
\mdensityCoeff_\ipart &:= \Prob[\mintervalBound_{\ipart - 1} \le \cintensity(\intime, \state) \le \mintervalBound_{\ipart}]\\ &= \int \int_0^\infty \Id_{(\mintervalBound_{\ipart - 1}, \mintervalBound_{\ipart}]}(\cintensity(\intime, \state)) \Density(\intime, \state) \idensity(\state) \dd \intime \dd \state.
\end{align*}

If we define $\intime^{(\cintensity)}_\itime(\statefrom)$ similarly as $\intime_\itime(\statefrom)$, namely by $\intime^{(\cintensity)}_0(\statefrom) = 0$ and recursively $\intime^{(\cintensity)}_\itime(\statefrom) := \min\{ \intime > \intime^{(\cintensity)}_{\itime - 1}(\statefrom): \cintensity(\intime, \statefrom) = \mintervalBound_\ipart \text{ for some } \ipart \}$, then

\begin{equation}
\mdensityCoeff_\ipart \approx \weight \cdot \sum_{\iPart = 1}^\inum \idensity(\statefrom_\iPart) \sum_{\itime: \iTrap(\itime,\statefrom_\iPart) =\ipart } \int_{\intime^{(\cintensity)}_{\itime-1}(\statefrom_\iPart)}^{\intime^{(\cintensity)}_{\itime}(\statefrom_\iPart)} \Density(\intime, \statefrom_\iPart) \dd \intime,
\label{eq:mdensityCoeffAppro}
\end{equation}

where $\iTrap(\itime, \statefrom) :=  \ipart$ if $\target(\intime, \statefrom) \in \statePartition_\ipart$ for $\intime_{\itime-1}(\statefrom) < \intime < \intime_\itime(\statefrom)$.
We further have

\begin{equation}
\Delta \cintensity \cdot \acid(\aci) \approx \sum_{\ipart = 1}^\infty \mdensityCoeff_\ipart\Id_{(\mintervalBound_{\ipart - 1}, \mintervalBound_{\ipart }]}(\aci).
\end{equation}

§ CASE STUDIES

§.§ Random telegraph permits an analytic solution

Continuing example <ref>, we consider the normalized parameters $\RTON := \RTOn/\gain, \RTOFF := \RTOff/\gain$ and, for convenience, drop the tilde again. For simplicity, $\cprob := \cprob_1$. Let $\RTRoot_1< \RTRoot_2$ be the roots of the quadratic equation at equilibrium

\begin{equation}
0 = \gain(\RTOn - (\RTOff + \RTOn + 1) \RTRoot + \RTRoot^2)
\label{eq:logistic_equilibrium}.
\end{equation}

Compute the solution of eq. (<ref>) that yields

\begin{equation}
\cprob(\intime) =  \RTRoot_2 - \frac{\Delta \RTRoot}{1 + \frac{1 - \RTRoot_1}{ \RTRoot_2 - 1} e^{-\gain\Delta \RTRoot \intime}}.
\label{eq:rt:cintensity}
\end{equation}

Then the unnormalized density $\Density(\intime)$ is obtained from solving eq. (<ref>)

\begin{equation}
\Density(\intime) =  e^{-\gain \RTRoot_1 \intime}\frac{\RTRoot_2 - 1}{\Delta\RTRoot} + e^{-\gain \RTRoot_2 \intime}\frac{1 - \RTRoot_1}{\Delta\RTRoot}
\label{eq:mixture_exponentials}
\end{equation}

and $\idensity = \E[\intensity_\infty] = \frac{\gain \RTOn}{\RTOn + \RTOff}$.

Using $\phi(\gain x) = \gain \phi(x) + x \phi(\gain)$ the mutual information rate can be written as

\begin{align}
\IR (\msg, \out) =  &- \gain \idensity  \int_0^\infty \phi(\cprob(\intime)) \Density(\intime) \dd \intime \label{eq:reparametrization}\\
= &- \frac{\gain\RTOn}{\RTOn + \RTOff} \int_0^\infty \phi\left(\RTRoot_2 - \frac{\Delta \RTRoot}{1 + \frac{1 - \RTRoot_1}{ \RTRoot_2 - 1} e^{-\Delta \RTRoot \intime}}\right) \times \nonumber\\
&\times \left[e^{- \RTRoot_1 \intime}\frac{\RTRoot_2 - 1}{\Delta\RTRoot} + e^{- \RTRoot_2 \intime}\frac{1 - \RTRoot_1}{\Delta\RTRoot} \right] \dd \intime.
\end{align}

The linear time scaling $\intime \mapsto \gain \intime$ was used in the second equality.
We identify eq. (<ref>) as a reparametrization of the integral representation previously reported [49] for $\RTRoot = \RTRoot_1$ and $\CProb_\ti(1) \overset{\mathrm{d.}}{\to} \condmean$

\begin{equation}
\IR (\msg, \out) = -\E[\gain\phi(Z)] = -\gain \int_{\RTRoot}^1 \phi(z) \pi_Z(z) \dd z.
\label{eq:integral_ISIT}
\end{equation}

The link can be obtained via the transformation rule applied on the transformation $\cprob\colon [0,\infty)\to (\RTRoot, 1]$, see Appendix <ref>.
The reparametrization (<ref>) has the advantage that vertical asymptotes of the integrand at $z = \RTRoot$ can be avoided. Furthermore, the integral bounds do not depend on the system parameters $\RTOn, \RTOff, \gain$, allowing more uniformly chosen integral bounds in the numerical approximation, compare remark <ref>.

§.§.§ Computation of partial derivatives

Another advantage lies in enabling the computation of partial derivatives. Let us pick up the tilde again. In order to answer questions of optimality of $\IR(\msg, \out)$, the partial derivatives of $\IR(\msg, \out) = \IR(\RTON, \RTOFF, \gain)$ are relevant. Note the relation $\partial_1 \IR(\RTON, \RTOFF, \gain) := \partial_{\RTON} \IR(\msg, \out) = \gain \partial_{\RTOn} \IR(\msg, \out)$. In particular the nullclines $[\partial_1 \IR(\RTON, \RTOFF, \gain) = 0]$ and $[\partial_2 \IR(\RTON, \RTOFF, \gain) = 0]$ contain information about optimal points. By the implicit function theorem, the $\RTON$-nullcline $\RTON \mapsto \graphNullcline(\RTON)$ satisfies the ODE

\begin{equation}
\graphNullcline_\gain'(\RTON) = -\frac{\partial_{11}\IR(\RTON, \graphNullcline(\RTON), \gain)}{\partial_{12}\IR(\RTON, \graphNullcline(\RTON), \gain)}.
\label{eq:IFT}
\end{equation}

Furthermore, to decide on convexity of the $\RTOn$-nullcline, positivity of
\begin{equation}
\graphNullcline_\gain''(\RTON) = \left[ \frac{2 \partial_{11}\IR \partial_{112}\IR}{\partial_{12}\IR^2} -\frac{\partial_{111}\IR}{\partial_{12}\IR} - \frac{\partial_{11}\IR^2 \partial_{122}\IR}{\partial_{12}\IR^3}\right](\RTON, \graphNullcline_\gain(\RTON), \gain)
\label{eq:convexity}
\end{equation}

must be checked.
This motivates to compute partial derivatives up to the third order.
The numerical method is exemplified for $\partial_1 \IR(\RTON, \RTOFF, \gain)$. In order to appreciate the reparametrization, we observe that the Leibniz rule for differentiation of the parameter integral in eq. (<ref>) fails when there is an asymptote at $\RTRoot_1$. Then the lower boundary term evaluates to $-\infty$. Thus, let us exploit the reparametrization.
Define $\cprob_1(\intime) := \partial_{\RTON} \cprob(\intime)$ and $\density_1(\intime) := \partial_{\RTON} \density(\intime)$ as well as

\begin{align*}
\cummIR_1(\intime) := \partial_{\RTON} \cummIR(\intime) &= \partial_{\RTON}  \int_0^\intime - \gain \phi(\cprob(\ti)) \density(\ti) \dd \ti\\
&=\partial_{\RTON}  \int_0^{\gain \intime} - \phi(\cprob(\ti/\gain)) \density(\ti/\gain) \dd \ti.
\end{align*}

For the evolution of $\cprob_1(\intime)$ one takes advantage of $\dot \cprob_1(\intime) = \partial_1\dot \cprob(\intime)$. The joint evolution of $\density, \cprob, \density_1, \cprob_1, \cummIR_1$ and initial values are given in the Appendix, <ref>.
The saturation value $\lim_{\intime \to \infty}\cummIR_1(\intime)$ is the partial derivative $\partial_1 \IR(\RTON, \RTOFF, \gain)$.

§.§.§ Mutual information in the phase plane

Phase plane analysis of the mutual information for the Poisson channel with random telegraph input. Nullclines $\partial_1 \bar{\I} = 0$ and $\partial_2 \bar{\I} = 0$ were evaluated, using $T = 1000$. Colored arrows indicate the gradient of the information rate, calculated alike. Optimization domains are rectangular. Depending on the location of the domain's upper right corner $(r_1, r_2)$, the optimum is assumed on the nullclines $\partial_1 \bar{\I} = 0$ [$(r_1, r_2)$ in C] and $\partial_2 \bar{\I} = 0$ [$(r_1, r_2)$ in A], respectively, or in the interior of region B [$(r_1, r_2)$ in B]. Regions I and II contain all $(r_1, r_2)$, whose optima favor the On and Off state, respectively. In the inset the dashed line $\RTOff = (e-1)\RTOn$ separates the nullclines. The gain was fixed at $\gain = 1$.

For the analysis we set $\gain = 1$, which is justified by the scaling behavior $\IR(\RTOn, \RTOff, \gain) = \gain \IR(\RTON, \RTOFF, 1) =: \gain \IR(\RTON, \RTOFF)$. For convenience we drop the tilde.
Along the line $(\RTOn, (e - 1)\RTOn)$, both partial derivatives are numerically found not to switch sign. From this we conclude that the nullclines do not intersect. This excludes local optima of $\IR(\RTOn, \RTOff)$. For a rectangular constraint $0 < \RTOn \le \OnConstraint, 0 < \RTOff \le \OffConstraint$, the maximizing pair $(\MaxOn, \MaxOff)$ is consequently always located on one boundary $\RTOn = \OnConstraint$ or $\RTOff = \OffConstraint$. The $(\RTOn, \RTOff)$-plane is split into regions $A, B$ and $C$ by the nullclines. The pair $[\sgn(\partial_1\I), \sgn(\partial_2\I)]$ characterizes the regions: $[1,-1],[1,1],[-1,1]$ on A,B,C. Depending on the location of $(\OnConstraint, \OffConstraint)$, the maximum $(\RTOn^\ast, \RTOff^\ast)$ shows different behavior. If located in A, it holds $\RTOn^\ast = \OnConstraint, \RTOff^\ast < \OffConstraint$, while a location in B enforces $\RTOn^\ast = \OnConstraint, \RTOff^\ast = \OffConstraint$, and finally, $\RTOn^\ast < \OnConstraint, \RTOff^\ast = \OffConstraint$ for $(\OnConstraint, \OffConstraint)$ in region C. As a summary, the maximum is always located in region B or its boundary. Numerically, we find that the $\RTOn$-nullcline crosses the bisection line $(\RTOn, \RTOn)$ at $\intersection$. An evaluation of $\partial_1 \IR(\RTOn, \RTOff)$ for $\RTOn \in (0, \intersection)$ indicates no further intersections. The derivative (<ref>) at $\RTOn = \intersection$ is found to be larger than one. Furthermore, eq. (<ref>) evaluates to positive values in the region $\{(\RTOn, \RTOff) \in [0, 0.3]^2: \RTOn - 0.07 \le\RTOff \le \RTOn + 0.01\}$. Hence, the isoclines that transit the region are convex, in particular the $\RTOn$-nullcline. The constraint $0 < \RTOff \le \OffConstraint$ for $\OffConstraint \in (0, \intersection)$ then returns a maximum $(\MaxOn, \MaxOff)$ with $\MaxOn > \MaxOff$. This shows that a bandwidth-like constraint can impose an On-favoring maximum complementing the classical result by Kabanov. The sequence of input processes for the Poisson channel that exhaust its capacity is a random telegraph process with $\RTOn, \RTOff \to \infty$ and $\RTOff/ \RTOn \to e-1$. The stationary input distribution favors the Off state, occupying it $1-1/e$ of the time.

§.§.§ Mutual information rate increases with channel gain

To answer the question whether $\partial_\gain\IR(\msg, \out)$ is positive, we first observe, see Appendix <ref>,

\begin{equation}
\partial_3\IR(\RTOn, \RTOff, \gain) =\partial_3\IR\left(\frac{\RTOn}{\gain}, \frac{\RTOff}{\gain}, 1\right).
\label{eq:scaling}
\end{equation}

Consequently, the analysis of $\partial_3\IR\left(\RTOn, \RTOff, 1\right)$ suffices.

\begin{align}
&\partial_3\IR\left(\RTOn, \RTOff, \gain\right)\vert_{\gain = 1} \nonumber \\ &= -\partial_\gain\left(\gain \E[\phi(\condmean)]\right)\vert_{\gain = 1}\nonumber\\
&= \IR\left(\RTOn, \RTOff, 1\right) -\partial_\gain\E[\phi(\condmean)]\vert_{\gain = 1}\nonumber\\
&= \IR\left(\RTOn, \RTOff, 1\right) - \partial_\gain \int_0^\infty \phi(\cprob(\ti)) \density(\ti) \dd \ti\vert_{\gain = 1}.
\label{eq:affirm}
\end{align}

The second term is evaluated analogously to ODE (<ref>) - (<ref>).
The numerical evaluation shows that the derivative with respect to the channel gain is non-negative.
We affirm a result in <cit.> using eq. (<ref>): For $\RTOn, \RTOff \to \infty$ with $\RTOn/\RTOff$ constant, $\condmean$ converges in distribution to a delta distribution at $\frac{\RTOn}{\RTOn + \RTOff}$ independent of $\gain$ and the second term of eq. (<ref>) vanishes. Consequently, we derive, using in the second equality eq. (<ref>) and (<ref>),

\begin{align*}
\partial_3\IR(\RTOn, \RTOff, \gain)\vert_{\gain = 0} &= \lim_{\gain \to 0} \partial_3\IR(\frac{\RTOn}{\gain}, \frac{\RTOff}{\gain}, 1) \\
&= \lim_{\RTOff \to \infty}\IR(\frac{\RTOn}{\RTOff}\RTOff, \RTOff, 1)\\
&= \E[\phi(\msg_\infty)] -\phi(\frac{\RTOn}{\RTOn+\RTOff})\\
&= \E[\phi(\msg_\infty)] - \phi(\E[\msg_\infty]).
\end{align*}

§.§ Random telegraph with dark current

Random telegraph input with dark current. ACIDs obtained from Monte Carlo samples and simulation-free computation agree. The density has a vertical asymptote at $\RTRoot_1$ and is non-smooth at $\jtarget_\infty = \jtarget(\RTRoot_1)$. Parameters were $\LeakageOn =  \LeakageOff = 0.1, \leakage = 0.1, \dynamicRange = 1$. Sample size was $10^6$ for the Monte Carlo simulation. Grid point number was $2000$ for the $\idensity$-support $(\jtarget_\infty, \ONState]$ and number of fixed point iterations was $\niter = 12$.

In the previous example, the dimension $\num$ of the sufficient statistic $\state$ was $0$ and no integral boundary condition was needed. We will proceed with a non-trivial $\state$ when continuing example <ref>. Similar as in <ref>, define $\tilde{\RTRoot}_1 < \tilde{\RTRoot}_2 $ as the roots of the quadratic equilibrium equation

\begin{equation}
0 = \RTOn \dynamicRange  - (\RTOff + \RTOn + \dynamicRange) \RTRoot + \RTRoot^2
\label{eq:leakage:logistic_equilibrium}
\end{equation}

and $\RTRoot_i :=\tilde{\RTRoot}_i + \leakage, i = 1,2 $. Obtain

\begin{equation}
\cintensity(\intime, \state) =  \RTRoot_2 - \frac{\Delta \RTRoot}{1 + \frac{\state - \RTRoot_1}{ \RTRoot_2 - \state} e^{-\Delta \RTRoot \intime}}
\label{eq:rt:cintensity}
\end{equation}

\begin{equation}
\Density(\intime, \state) =  e^{- \RTRoot_1 \intime}\frac{\RTRoot_2 - \state}{\Delta\RTRoot} + e^{- \RTRoot_2 \intime}\frac{\state - \RTRoot_1}{\Delta\RTRoot}.
\label{eq:RT_survival_time}
\end{equation}

The number of solutions of $\target(\intime, \statefrom) = \stateto$ is $N(\stateto, \statefrom) \in \{0,1\}$ with

\begin{equation}
\intime(\stateto, \statefrom) = \Delta \RTRoot^{-1}\left\{\ln\left(\frac{\RTRoot_2 - \jtarget^{-1}(\stateto)}{\jtarget^{-1}(\stateto) - \RTRoot_1}\right) + \ln\left(\frac{\statefrom - \RTRoot_1}{\RTRoot_2 - \statefrom}\right)\right\}.
\end{equation}

Observe that for every $\state$ the trajectory $\cintensity(\cdot, \state)$ is decreasing and $\lim_{\intime \to \infty} \cintensity(\intime, \state) = \RTRoot_1$. Consequently, for any $\intime, \state$

\begin{equation}
\target(\intime, \state) > \jtarget(\RTRoot_1) = \frac{(\RTRoot_1 - \leakage)(\ONState - \RTRoot_1)}{\RTRoot_1} =: \jtarget_\infty.
\end{equation}
We partition the $\idensity(\state)$-support $\stateSupport = (\jtarget_\infty, \ONState]$ into equidistant intervals $(\intervalBound_{\ipart-1}, \intervalBound_{\ipart}]$ with $\intervalBound_{\ipart} = \jtarget_\infty + \ipart \cdot \frac{1 - \jtarget_\infty}{\inum}, \ipart = 1, \dots, \inum$. Choose representatives $\state_\ipart = \frac{\intervalBound_{\ipart} + \intervalBound_{\ipart-1}}{2}$. Then the matrix entries in (<ref>) are given by

\begin{equation}
\OpI_{\ipart, \iPart} = \Density(\intime(\intervalBound_{\ipart - 1}, \state_\iPart)\vee 0, \state_\iPart) - \Density(\intime(\intervalBound_{\ipart} , \state_\iPart) \vee 0, \state_\iPart)
\end{equation}

\begin{equation}
\Density(\intime(\stateto, \statefrom), \statefrom) = \left( \frac{\jtarget^{-1}(\stateto) - \RTRoot_1}{\statefrom - \RTRoot_1} \right)^{\frac{\RTRoot_1}{\Delta \RTRoot}}\left( \frac{\RTRoot_2 - \statefrom}{\RTRoot_2 - \jtarget^{-1}(\stateto)} \right)^{\frac{\RTRoot_2}{\Delta \RTRoot}}.
\end{equation}

The initial condition $\idensity(\state)$ was found by fixed-point iteration $\OpI^{2^\niter}$ of (<ref>) with $\niter$ iterations. The mutual information rate was computed via (<ref>) and (<ref>).
For (<ref>) we compute

\begin{equation}
\intime^{(\cintensity)}(\stateto, \statefrom) = \Delta \RTRoot^{-1} \left\{ \ln\left(\frac{\RTRoot_2 - \stateto}{\stateto - \RTRoot_1}\right) - \ln\left(\frac{\RTRoot_2 - \statefrom}{\statefrom - \RTRoot_1}\right) \right\}
\end{equation}

and explicitly obtain in (<ref>)

\begin{align}
&\int_{\intime^{(\cintensity)}(\mintervalBound_{\ipart}, \statefrom)}^{\intime^{(\cintensity)}(\mintervalBound_{\ipart-1}, \statefrom)} \Density(\intime, \statefrom) \dd \intime =  \frac{(\RTRoot_2 - \statefrom)^{\frac{\RTRoot_2}{\Delta \RTRoot}}}{ \Delta \RTRoot(\statefrom - \RTRoot_1)^{\frac{\RTRoot_1}{\Delta \RTRoot}}}\times\nonumber\\
&\qquad \qquad \left[ \frac{(\cintensity - \RTRoot_1)^\frac{\RTRoot_1}{\Delta \RTRoot}}{\RTRoot_1(\RTRoot_2 - \cintensity)^\frac{\RTRoot_1}{\Delta \RTRoot}} + \frac{(\cintensity - \RTRoot_1)^\frac{\RTRoot_2}{\Delta \RTRoot}}{\RTRoot_2(\RTRoot_2 - \cintensity)^\frac{\RTRoot_2}{\Delta \RTRoot}}\right]_{\cintensity = \mintervalBound_{\ipart - 1}}^{\cintensity = \mintervalBound_{\ipart}} .
\end{align}

Fig. <ref> shows the agreement of the simulation-free computation with the Monte-Carlo sampling. The ACID has an asymptote at $\RTRoot_1$ and is not differentiable at $\jtarget_\infty = \jtarget(\RTRoot_1)$.

Results for different $\RTOn, \RTOff$ and increasing dark current are shown in figure <ref>. The information rate was computed with ODE system (<ref>) - (<ref>). For any examined pair $(\LeakageOn, \LeakageOff)$, the information rate decreases with dark current as expected. The figure reveals a notable property. For fixed $\LeakageOn = 0.1$, the plots for $\LeakageOff = 0.1$ and $\LeakageOff = 0.5$ intersect. This means that increasing the dark current increases the information rate $\IR(0.1, 0.1)$ relative to $\IR(0.1, 0.5)$. Consequently, an increased dark current can qualitatively alter monotonicity and optimality properties in the $(\LeakageOn, \LeakageOff)$-phase plane. For example, the On favoring region increases with dark current [49].

Random telegraph input with dark current $\leakage$ and amplitude $\dynamicRange = 1$. The information rate $\IR(\RTOn, \RTOff)$ is plotted for different values of $\leakage, \RTOn, \RTOff$. The graphs of $\RTOn = 0.1, \RTOff = 0.1$ and $\RTOn = 0.1, \RTOff = 0.5$ intersect, showing that dark current can alter the monotonicity properties of $\IR(\RTOn, \RTOff)$ in the $( \RTOn,\RTOff)$-plane.

§.§ Double On single Off

Binary Markov input processes exhaust the capacity when their switching rates tend to infinity. The defining property for optimality in the limit is only the proportion On/Off. With the autocorrelation time going to $0$ for an exhausting sequence of Semi-Markov processes also, the Markov property might as well be relaxed.
Consider the following capacity problem: We restrict the input process class to binary semi-Markov processes and impose a lower bound on the average sojourn times in the On and the Off state. Is the Markov case with its exponential sojourn times the capacity-achieving input? Here, we consider a Semi-Markov processes with exponential sojourn time in the Off and Erlang sojourn time in the On state. This can be realized as an instance of example class <ref>, the afore-mentioned three-state Markov process that circles through one inactive and two active states.
For the numerical evaluation we discretize $[0,1] \ni \state$. The state variables $( \Donsofffrac(\intime, \state), \z(\intime, \state))$ evolve according to (<ref>) - (<ref>)
with initial conditions $[ \Donsofffrac(0, \state), \z(0, \state)] = [\state, 1]$ and $\target(\intime, \state) = \Donsofffrac(\tau, \state)$. The times $\intime_\itime(\statefrom_\iPart)$ in (<ref>), that satisfy

\begin{equation}
\target(\intime, \statefrom_\iPart) = \intervalBound_{\ipart},
\label{eq:Donsoff_event}
\end{equation}

were found by evolving the ODE system (<ref>), (<ref>), (<ref>) and checking for the event (<ref>). The matrix entries (<ref>) were evaluated and $\idensity(\state)$ was found by fixed-point iteration with $2^\niter$ iterations. The rates $\DonsoffOn = \OnValue, \DonsoffR = \DonsoffOff = \OffValue$ exhibit a mutual information of $\IR(\msg, \out) = \DonsoffMIValue$ compared to $\IR(\msg, \out) = \RTMIValue$ for the Markov case.
The rates were tuned such that the average sojourn times are the same. This example shows that the Markov input does not generally solve the capacity problem with average sojourn time constraint. It remains an open research question which On and Off sojourn time distributions are capacity-achieving.

The relevance for non-Markovian binary models is also reflected in biological systems. For instance, promoter models can exhibit multiple active and inactive states [50] while $\out_\ti$ counts mRNA synthesis events.

The Double On Single Off model increases the mutual information rate compared to the random telegraph model. Parameters were $\DonsoffOn = \RTOn = \OnValue$ and $\DonsoffR = \DonsoffOff = 2\cdot \RTOff$. This choice guaranteed a model match in terms of mean On and Off times.

§.§ Hawkes process

The Hawkes process, introduced in <ref>, is fully characterized in terms of second-order properties [19], which was used in [8] to obtain capacity upper bounds via the link to optimal linear estimation. Beyond information theoretic applications, the usefulness of the Hawkes process' ACID has been mentioned in the literature [20, 17]. Oakes [20] suggested that its equilibrium distribution, i.e., the ACID, can be found by iteratively solving an integral equation, but does not elaborate more on it. Daley & Vere-Jones <cit.> provides a manual for setting up the integral equation and solving it with the method of steps. We are not sure whether Oakes meant this method by "iterative solution". We provide a different approach, using BReT-P and the fixed point iteration.
The support of $\idensity(\state)$ is $(\HawkesBase + \HawkesJump, \infty)$ and contained in $(\HawkesBase, \infty)$. We consider the equidistant partition $\intervalBound_\ipart = \HawkesBase  + \ipart \Delta\state, \ipart = 0, \dots, \inum$ for $\Delta \state\cdot \inum$ large enough to cover most of the probability weight, i.e., for

\begin{equation}
\int_0^{\Delta \state\cdot \inum} \idensity(\state + \HawkesBase + \HawkesJump) \dd \state \approx \int_0^{\infty} \idensity(\state + \HawkesBase + \HawkesJump) \dd \state.
\end{equation}

\begin{equation}
\intime(\stateto, \statefrom) = \HawkesDecay^{-1}\ln\left( \frac{\statefrom - \HawkesBase}{\stateto - \HawkesBase - \HawkesJump} \right)
\end{equation}

\begin{equation}
\Density(\intime(\stateto, \statefrom), \statefrom) = \left(\frac{\stateto - \HawkesJump - \HawkesBase}{\statefrom - \HawkesBase} \right)^{\frac{\HawkesBase}{\HawkesDecay}} e^{-\HawkesDecay^{-1}(\statefrom - \stateto + \HawkesJump)}.
\end{equation}

Examples of the ACID, obtained from fixed point iteration as described in <ref>, are shown in fig. <ref>. These may serve as initial distribution of $\CIntensity_0$ for the stationary Hawkes process. Using martingale theory, the equilibrium variance of $\CIntensity_\ti$ was derived (see Appendix <ref>) to equal

\begin{equation}
\V[\CIntensity_\infty] = \frac{\HawkesDecay\HawkesBase\HawkesJump^2}{2(\HawkesDecay - \HawkesJump)^2}.
\label{eq:HawkesVariance}
\end{equation}

The parameter sets in fig. <ref> were chosen, such that the ACID's first and second order moments are constant, but vary in the exponential decay parameter $\HawkesDecay$ of the Hawkes kernel. This makes the shape vary qualitatively. While for fast decay ($\HawkesDecay$ large) the region near the base value $\HawkesBase$ is frequented more heavily, for slow decay, the CI spends more time in the middle regime around the mean $\gain\HawkesMean$. This illustrates that the ACID analysis goes beyond the mean and variance analysis, i.e., that the ACID is parameterized by more than two parameters.

§.§.§ The Hawkes approximate marginal simulation

For the purpose of comparing the Hawkes process to Markov-modulated Poisson processes and approximate filters, we view it as an approximate marginal simulation. It is obtained from the optimal linear filter approximation of the following process class. Let $\out_\ti$ be a doubly stochastic Poisson process whose external signal $\msg_\ti$ has the following first- and second-order statistics

\begin{equation}
\E[\msg_\ti] \equiv \HawkesMean, \quad \Cov[\msg_\ti, \msg_s] = \HawkesVariance e^{-\HawkesAutocov \vert \ti - s \vert}.
\label{eq:HawkesFirstSecond}
\end{equation}

And the $\sigmalg_\ti^{\msg, \out}$-intensity is $\intensityFct_\ti = \gain \msg_\ti$. External signals $\msg_\ti$ of this form comprise
(i) the random telegraph model with or without dark current and
(ii) the CIR process.
The optimal linear filter theory by Snyder [3] identifies the estimator $\CIntensity_\ti$ that minimizes $(\tilde\intensity_\ti - \intensity_\ti)^2$ among all estimators of the following form, which is linear in $\out_{[0, \ti]}$,

$$\tilde\intensity_\ti = a(\ti) + \int_0^\ti h(\ti, u) \dd \out_u. $$

The resulting form of $\CIntensity_\ti$ appeals as a variant of the Kalman filter, with the Riccati equation for the innovation gain:

\begin{align}
\dd \CIntensity_\ti = -\HawkesAutocov(\CIntensity_\ti - \HawkesMean) \dd \ti + \HawkesJump(\ti) (\dd \out_\ti - \CIntensity_\ti \dd \ti)\label{eq:Kalman}\\
\frac{\dd}{\dd \ti} \HawkesJump(\ti) = -\HawkesJump(\ti)^2 - 2\HawkesAutocov \HawkesJump(\ti) + \frac{2\gain\HawkesAutocov\HawkesVariance}{\HawkesMean}. \nonumber
\end{align}

At equilibrium $\HawkesJump(\ti)$ can be replaced by the constant $\HawkesJump$ that solves $\HawkesJump^2 + 2\HawkesAutocov \HawkesJump - \frac{2\gain\HawkesAutocov\HawkesVariance}{\HawkesMean} = 0$. We interpret the evolution equation (<ref>) as described in section <ref>. For the approximate marginal simulation $\out_\ti$ is self-exciting with CI $\CIntensity_\ti$ in contrast to $c\E[\msg_\ti\vert \sigmalg^\out_\ti]$ for estimation. So the eq. (<ref>) is of the same shape as the Hawkes process.
A parameter match links the quadruple $(\HawkesMean, \HawkesVariance, \HawkesAutocov, \gain)$ to the original triple $(\HawkesBase, \HawkesJump, \HawkesDecay)$ as follows

\begin{equation}
\HawkesJump = \sqrt{\HawkesAutocov^2 + \frac{2 \gain \HawkesAutocov\HawkesVariance}{\HawkesMean}} - \HawkesAutocov, \; \HawkesDecay = \HawkesAutocov + \HawkesJump,\; \HawkesDecay\HawkesBase = \gain\HawkesMean\HawkesAutocov.
\label{eq:HawkesParameterLink}
\end{equation}

In the literature eq. (<ref>) are also called the mean intensity and covariance density of the Hawkes process [51].

We emphasize that the stationary distributions for common input processes, such as the CIR process, birth-death process or random telegraph model, are entirely characterized by mean and variance. The decay parameter $\HawkesAutocov$ in eq. (<ref>) is not captured by the stationary distribution. In contrast, the ACID does capture a change in $\HawkesAutocov$, see fig. <ref>, so it contains temporal information about the input process.

title = $\HawkesDecay = 0.3$,
name = plot_left,

scale only axis,
xlabel = $\CIntensity_\infty$,
axis background/.style=fill=white
[color=mycolor1, line width = 1pt]
table[row sep=crcr]
0.396706523493108	0.0289978038965

0.411717553273543	0.0443344905540197

0.426728583053978	0.0543757833620839

0.441739612834413	0.0632155471606231

0.501783731956151	0.0968999530184389

0.516794761736587	0.105736088444382

0.531805791517021	0.114900114762284

0.546816821297456	0.124444958045024

0.561827851077891	0.134417738428697

0.576838880858325	0.144862217926052

0.59184991063876	0.155820356173325

0.606860940419194	0.167333352909521

0.62187197019963	0.179442381460333

0.636882999980064	0.191154741057075

0.651894029760498	0.201720409751359

0.666905059540934	0.211741850197201

0.696927119101803	0.231303833899674

0.726949178662673	0.250541182123016

0.756971238223542	0.269504835027487

0.786993297784411	0.288094769040295

0.802004327564847	0.297192561134054

0.817015357345281	0.30611852335047

0.832026387125715	0.314836382149893

0.84703741690615	0.323305010256004

0.862048446686585	0.331478271397654

0.87705947646702	0.339332537374218

0.892070506247454	0.346902526318541

0.907081536027889	0.35423778105089

0.922092565808324	0.36135597537213

0.937103595588758	0.368258895633782

0.952114625369194	0.374944108545156

0.967125655149628	0.381407218841238

0.982136684930063	0.387642941487885

0.997147714710498	0.393645766453441

1.01215874449093	0.399410449503081

1.02716977427137	0.404932423166543

1.0421808040518	0.410208173209008

1.05719183383224	0.415235605489414

1.07220286361267	0.420014418469905

1.08721389339311	0.424546491740234

1.10222492317354	0.428836298293247

1.11723595295398	0.432890661610737

1.13224698273441	0.436716600844704

1.14725801251485	0.440318830522513

1.16226904229528	0.4436999052044

1.17728007207571	0.446861695636224

1.19229110185615	0.449806052951572

1.20730213163658	0.452535017699837

1.22231316141702	0.455050928882234

1.23732419119745	0.457356482316603

1.25233522097789	0.459454755431433

1.26734625075832	0.461349205420258

1.28235728053876	0.463043643399656

1.29736831031919	0.464542185039413

1.31237934009963	0.465849176851416

1.32739036988006	0.466969096461134

1.3424013996605	0.467906424528368

1.35741242944093	0.468665501629443

1.37242345922137	0.469250437757964

1.3874344890018	0.469665152717921

1.40244551878224	0.469913503250726

1.41745654856267	0.469999375948633

1.4324675783431	0.469926716511656

1.44747860812354	0.469699531873804

1.46248963790397	0.46932188237638

1.49251169746484	0.4681316126746

1.52253375702571	0.466388872449148

1.55255581658658	0.464126399234116

1.58257787614745	0.461376189488973

1.61259993570832	0.458169244048023

1.64262199526919	0.454535688685344

1.67264405483006	0.450504999883778

1.70266611439093	0.446106002768364

1.7326881739518	0.441366742916623

1.76271023351267	0.436314363698751

1.79273229307354	0.43097502830033

1.82275435263441	0.425373895671092

1.85277641219528	0.419535138478663

1.88279847175615	0.413481969983473

1.92783156109745	0.404048621070174

1.97286465043875	0.394253728240344

2.01789773978006	0.384163392151796

2.0779418589018	0.37035523551494

2.15299700780397	0.352709055490348

2.4231955438518	0.288693438885065

2.48323966297353	0.274868771394396

2.54328378209527	0.261322002836589

2.60332790121701	0.248096875762265

2.66337202033875	0.235229339985134

2.70840510968006	0.225830998968971

2.75343819902136	0.21666004226056

2.79847128836266	0.207724480335622

2.84350437770397	0.199030646948521

2.88853746704527	0.190583345700784

2.93357055638658	0.182385988518453

2.97860364572788	0.174440726073072

3.02363673506919	0.166748570275426

3.06866982441049	0.159309509025645

3.11370291375179	0.152122613430826

3.1587360030931	0.145186137726749

3.2037690924344	0.138497612168948

3.24880218177571	0.132053929179198

3.29383527111701	0.12585142304494

3.33886836045831	0.119885943475128

3.38390144979962	0.114152923320032

3.42893453914092	0.10864744076372

3.47396762848223	0.103364276296086

3.53401174760397	0.0966564098061413

3.5940558667257	0.0903205542509955

3.65409998584744	0.0843427776891899

3.71414410496918	0.078708855423276

3.77418822409092	0.0734044062836361

3.83423234321266	0.068415010480666

3.8942764623344	0.0637263107851282

3.95432058145614	0.0593240986828905

4.01436470057788	0.0551943870362006

4.07440881969962	0.0513234706714698

4.13445293882136	0.0476979762042982

4.20950808772353	0.0434914640317441

4.2845632366257	0.0396236396638754

4.35961838552788	0.0360710605116283

4.43467353443005	0.0328113883010177

4.50972868333222	0.0298234167318068

4.5847838322344	0.0270870822251128

4.67485001091701	0.0241090760923779

4.76491618959961	0.0214356692305895

4.85498236828222	0.0190385654340064

4.96005957674526	0.0165564490009693

5.06513678520831	0.0143770933307223

5.17021399367135	0.0124662752966431

5.29030223191483	0.0105719136728029

5.42540149993874	0.00876115814428324

5.56050076796265	0.00724072042308599

5.710611065767	0.00583848026224665

5.87573239335178	0.00458607708022996

6.055864750717	0.00350159196061472

6.26601916764309	0.00253003105059513

6.50619564413004	0.00171497537833787

6.7914052099583	0.00104384402267765

7.1366588949083	0.000524560740059243

7.58698978832134	0.00014436341762547

7.88721038393004	2.80399556018551e-06

title = $\HawkesDecay = 1$,
name = plot_middle,
at=(plot_left.right of south east), anchor=left of south west,
scale only axis,
xlabel = $\CIntensity_\infty$,
yticklabel = ,
axis background/.style=fill=white
[color=mycolor1, line width = 1pt, forget plot]
table[row sep=crcr]
0.772984443736583	0.88479860875823

0.791089286209328	0.6295823645754

0.809194128682073	0.566244126761044

0.82729897115482	0.532066594456809

0.845403813627565	0.510388696951587

0.86350865610031	0.495575702444857

0.881613498573055	0.485065604979575

0.8997183410458	0.477489859635451

0.917823183518545	0.472037467773648

0.935928025991291	0.468191397623723

0.954032868464036	0.465603596567689

0.972137710936781	0.464029683245432

0.990242553409526	0.463292211012778

1.00834739588227	0.463258765246351

1.02645223835502	0.463828274007795

1.04455708082776	0.464922111187267

1.06266192330051	0.46647812075752

1.08076676577325	0.468446488458863

1.098871608246	0.470786819470202

1.11697645071874	0.473466025274245

1.13508129319149	0.476456766715499

1.15318613566424	0.479736287581515

1.18939582060973	0.487088439047193

1.22560550555522	0.495403059920676

1.26181519050071	0.504594359543157

1.2980248754462	0.514600225177244

1.33423456039169	0.525375714434446

1.37044424533718	0.53688866557783

1.38854908780992	0.53406260180072

1.40665393028267	0.52409926440793

1.42475877275541	0.516637615912515

1.44286361522816	0.510438272371642

1.4609684577009	0.504968112975126

1.47907330017365	0.499988230724163

1.51528298511914	0.491000791365256

1.55149267006463	0.482838857195668

1.60580719748287	0.471432701591285

1.76875077973757	0.43809191689094

1.82306530715581	0.42613969244309

1.8592749921013	0.417722871955508

1.89548467704679	0.408873845249765

1.93169436199228	0.399532589992869

1.96790404693777	0.389640493544347

2.0584282593015	0.364001013153798

2.11274278671974	0.349558336816523

2.18516215661072	0.330966493260977

2.2575815265017	0.312895623141621

2.34810573886542	0.290883448012186

2.43862995122915	0.269500700819535

2.51104932112013	0.252930668542442

2.58346869101111	0.236970595852354

2.63778321842935	0.225501805099276

2.69209774584759	0.214489589721481

2.74641227326582	0.203896517988026

2.8188316431568	0.190389849269996

2.89125101304778	0.177567157791325

2.96367038293877	0.16541704803929

3.03608975282975	0.153933618834905

3.10850912272073	0.143110845142687

3.18092849261171	0.132937756370854

3.25334786250269	0.123393782621033

3.32576723239367	0.114447702344568

3.39818660228465	0.106068888117271

3.47060597217564	0.0982306909261013

3.54302534206662	0.0909079937523405

3.6154447119576	0.0840760082915608

3.68786408184858	0.0777097057439704

3.76028345173956	0.0717836789855841

3.85080766410329	0.0649564139394112

3.94133187646701	0.0587286658335504

4.03185608883074	0.0530549474944486

4.12238030119447	0.0478925653164843

4.21290451355819	0.0432010736163111

4.32153356839467	0.0381390856814505

4.43016262323114	0.033638247300253

4.53879167806761	0.0296416359990026

4.66552557537683	0.0255471361265904

4.79225947268604	0.0219936932976044

4.937098212468	0.0185096628013657

5.08193695224997	0.0155569745172812

5.24488053450468	0.0127750659267836

5.42592895923213	0.0102450937626646

5.62508222643233	0.00801995248346543

5.84234033610527	0.00612489217262002

6.0958081307237	0.00445784570117702

6.38548561028763	0.00308724787175407

6.72947761726979	0.00198296626366812

7.16399383661568	0.00112037337344617

7.76145363821627	0.000495748745336755

8.66669576185354	0.000124764953694623

9.80730083763649	4.0690877156635e-07

title = $\HawkesDecay = 3$,
name = plot_right,
at=(plot_middle.right of south east), anchor=left of south west,
scale only axis,
xlabel = $\CIntensity_\infty$,
yticklabel = ,
axis background/.style=fill=white
[color=mycolor1, line width = 1pt, forget plot]
table[row sep=crcr]
1.14331513983232	6.20152105104368

1.16697960313961	1.87039427343463

1.19064406644691	1.35425380446762

1.2143085297542	1.10322538794085

1.2379729930615	0.949564563563168

1.26163745636879	0.844065161199959

1.28530191967609	0.766384442534846

1.30896638298338	0.706414812390891

1.33263084629067	0.65850650772083

1.35629530959797	0.619229170253835

1.37995977290526	0.586367656433916

1.40362423621256	0.558421536969552

1.42728869951985	0.534334915315494

1.47461762613444	0.49486875254191

1.52194655274903	0.463844121576466

1.56927547936362	0.438796894926256

1.61660440597821	0.418156287503773

1.68759779590009	0.393245356116948

1.75859118582197	0.373625256790774

1.82958457574386	0.357848982289932

1.92424242897304	0.34119605728365

2.04256474550951	0.325461970538324

2.16088706204598	0.313760869320889

2.30287384188975	0.303529262316292

2.42119615842622	0.297371431791875

2.44486062173351	0.27560581153652

2.46852508504081	0.247430803904489

2.4921895483481	0.235392842339852

2.53951847496269	0.21840972121857

2.61051186488458	0.199916084278664

2.70516971811375	0.18103268693552

2.82349203465023	0.162139251154242

2.98914327780129	0.140464613139081

3.20212344756694	0.117115426015687

3.50976147056177	0.087958764160609

3.8173994935566	0.0623452159316038

4.05404412662954	0.0492425101903677

4.36168214962437	0.0362691664576857

4.74031356254108	0.0244573249218618

5.14260943876509	0.0159698196679017

5.68689209483286	0.00901353968925278

6.3968259940517	0.00421783713372648

7.48539130618724	0.00130171946520896

9.75717978368751	0.000105272201100348

12.9518823301723	3.8941706748119e-08

[anchor = north] at (plot_left.above north west) a);
[anchor = north] at (plot_middle.above north west) b);
[anchor = north] at (plot_right.above north west) c);

ACID for the Hawkes process. ACID mean $\gain\HawkesMean = \frac{\HawkesDecay\HawkesBase}{\HawkesDecay - \HawkesJump} = 2$ and ACID variance $\frac{\HawkesDecay\HawkesBase\HawkesJump^2}{2(\HawkesDecay - \HawkesJump)^2} = 1$ were constant, while $\HawkesDecay \in \{0.3, 1, 3 \}$ varied. The truncation $\Delta \state\cdot \inum$ of the support was chosen to be the $0.999$-quantile of the Gamma distribution with mean $\gain\HawkesMean$ and variance $\gain^2\HawkesVariance$. Discretization granularity was $\inum = 200$. Additionally, three equidistant representatives $\state_\ipart$, see eq. (<ref>) were chosen in each interval. Their mean function evaluations were used in eq. (<ref>) for the coefficients of $\OpI$. Number of iterations was $L = 15$.

We proceed with the second order analysis on the level of the output $\out_\ti$ in contrast to the input $\lambda_\ti$, and observe that for the process class in eq. (<ref>) the asymptotic mean and variance are known

$$\lim_{\ti \to \infty} \frac{1}{\ti} \E[\out_\ti] = \gain\HawkesMean, \quad \lim_{\ti \to \infty} \frac{1}{\ti} \V[\out_\ti] = \gain\HawkesMean + \frac{2\gain^2\HawkesVariance}{\HawkesAutocov}.$$

They are shared with the Hawkes process, see Appendix <ref>. Hence, the exact and the approximate marginal simulation cannot be discriminated by first and second order analysis of $\out_\ti$. ACID is employed to detect the difference.

§.§.§ Comparing the ACIDs of random telegraph with dark current and of the Hawkes process

The random telegraph with dark current modulates a Poisson process. Comparison of the exact Snyder filter and the optimal linear filter (Hawkes process). a), d) Lower panel shows 100 realizations of the approximate marginal simulation (blue, increasing) and the exact simulation (red, decreasing). The lower y axis was flipped to simplify visual comparison. Upper panel compares $\CIntensity_\ti$ for the Hawkes (blue) and the exact (red, dashed). The trajectory $\out_\ti$ for which both where computed, is shown in the lower panel (black). b), c) The cdfs of the ACIDs are compared for Hawkes (blue) and Snyder filter (red, dashed). e) The Wasserstein metric between the Hawkes ACID and the exact ACID is depicted for increasing dark current. The dashed line depicts the path metric (<ref>) obtained from a) as the average distance between both paths. (The trajectory $\out_\ti$ was obtained from exact simulation.) All plots used the parameters $\dynamicRange = 1, \LeakageOn = \LeakageOff = 0.1, \gain = 1$. Dark current varied in e), while a), b) used $\leakage = 0.1$,  and c), d) used $\leakage = 4$. Both values are indicated as dots in e).

For demonstration purposes we consider the tractable random telegraph input with dark current. It belongs to the considered input process class having exponentially decaying autocovariance function with $\HawkesAutocov = \LeakageOn + \LeakageOff$.
We compare the optimal linear filter (i.e. the Hawkes process) to the exact filter, obtained in <ref>. As discussed at the end of the previous paragraph, a first and second analysis cannot reveal a difference. Also a visual inspection of a sample trajectory ensemble can hardly tell them apart (lower panels in fig. <ref>a, d). And an asymptotic distribution for $\out_\ti$ cannot be compared because it does not exist. ACID is used to detect the parameter regimes where the approximate marginal simulation deviates from the exact marginal simulation (fig. <ref>b, c). Figure <ref>e) shows that for fixed switching rates, gain and amplitude, the deviation gets more severe for smaller dark current.

While on the one hand two counting processes that agree in the path distribution share the ACID, on the other hand the same or similar ACID does not necessarily imply a similar path measure. ACID only gives the ensemble picture of the CI value at a typical time point after the process entered stationarity. To have a more path-wise comparison, the following metric is considered:

\begin{equation}
\lim_{T \to \infty} \frac 1 T \int_0^T \vert\CIntensity_\ti - \CIntensity^H_\ti \vert \dd \ti .
\label{eq:pathwise_metric}
\end{equation}

Here, $\CIntensity_\ti$ is the exact filter and $\CIntensity^H_\ti$ is the optimal linear filter, both evaluated as functions of the history $\out_{[0,\ti]}$ of a drawn sample path $\out_{[0,\infty]}$.
Ergodicity guarantees that the metric does not depend on the sample path. However, it is not clear whether to simulate the sample path from the exact process with CI $\CIntensity_\ti$ or the approximate process with CI $\CIntensity^H_\ti$. (For a sample trajectory obtained from $\CIntensity_\ti$ it is shown in fig. <ref>e.) It depends on the context, i.e., the approximation goal, when to use a comparison of the ACIDs and when to use a path-wise comparison of the CIs for assessing an approximation. When approximating the information rate via eq. (<ref>) with an approximate filter $\CIntensity_\ti$, a comparison of the ACIDs seems suited.

Figure <ref> was obtained via our simulation-free method. The truncation $\Delta \state\cdot \inum$ of the support was chosen to be the $0.999$-quantile of the Gamma distribution with mean $\gain\HawkesMean$ and variance $\gain^2\HawkesVariance$. This choice was educated by the optimal linear approximation property and the fact that the stationary distribution of $\intensity_\ti$ is the above Gamma distribution. All subfigures in Fig. agree in the ACID variance (<ref>) and ACID mean. But the exponential decay of the Hawkes kernel varies. This makes the shape vary qualitatively. While for fast decay ($\HawkesDecay$ large), the region near the base value $\HawkesBase$ is frequented more heavily, for slow decay, the CI spends more time in the middle regime around the mean $\gain\HawkesMean$. This illustrates that the ACID analysis goes beyond the mean and variance analysis.

§.§ Gamma filter

The Gamma filter is used to illustrate the method for $\num = 2$. For technical details, we refer to Appendix <ref>.

§.§.§ Comparing the ACIDs of Gamma filter and Hawkes process

Consider the Markov-modulated Poisson process $\out_\ti$ whose external signal is a CIR-process, i.e., satisfies eq. $\eqref{eq:HawkesFirstSecond}$. The Hawkes process seen as optimal linear filter and the Gamma filter both approximate the true CI of $\out_\ti$. We compare them for the same parameters $\HawkesMean, \HawkesVariance, \HawkesAutocov, \gain$.
First, we inspect the asymptotic mean slope and variance slope for the Gamma filter. By martingale techniques (see Appendix <ref>), we can show that

\begin{equation}
\lim_{\ti \to \infty} \frac{1}{\ti}\V[\out_\ti^G] = \gain\GammaMu + \frac{2\gain^2(\E[\GammaVar_{\infty}] + \V[\GammaMean_\infty])}{\GammaDecay }.
\label{eq:FanoGamma}
\end{equation}

The Gamma filter satisfies the variance decomposition (see Appendix, A)

\begin{equation}
\E[ \GammaVar_{\infty}] + \V[ \GammaMean_\infty] = \GammaSigma,
\label{eq:VarianceDecomposition}
\end{equation}

consequently, the Gamma filter agrees with the exact CIR-modulated Poisson process in asymptotic first and second order moment. First and second order analysis cannot tell the Hawkes and Gamma filter apart.
A comparison of the ACIDs for a range of parameters reveals a slight increase of the Wasserstein metric for increased $\HawkesVariance$. The parameter $\HawkesAutocov$ had little effect on the Wasserstein metric. For technical details, see Appendix <ref>.
Fig. <ref>b (Appendix) depicts the example with the largest Wasserstein metric among the considered parameters, revealing that the ACIDs are still very similar. Due to the ACID's limitation as partial characteristic, we cannot deduce that the path measures are close in some notion of distance.
However, we conclude the following. When quantities are computed that only depend on the ACID, the optimal linear filter - appealing with efficient analytic expressions - might replace the Gamma filter.

The mutual information rate along the Poisson channel was efficiently approximated by Monte Carlo simulation in <cit.> via the Gamma filter. We replaced the Gamma filter by even the more efficient optimal linear filter.
A comparison for this case is shown in Appendix, <ref>.

Note, that both Gamma and Hawkes yield only an approximation of the exact information rate.
Between them, the Hawkes can be preferred in this case with its gain in efficiency and no loss in accuracy relative to the Gamma.
In what respect the replacement works for more complicated reaction networks, must be carefully evaluated.

§.§ Direct method for Hawkes and Dark Current

In special cases, the fixed point method can be applied directly to the ACID $\acid(\aci)$, i.e., we derived a linear fixed point equation

\begin{equation}
\acid(\aci) = \int \aciKernel(\aci, \acifrom)\acid(\acifrom) \dd \acifrom
\label{eq:aciKernel}
\end{equation}

for examples <ref> and <ref>.

Let $\PDMP_\ti$ follow the stochastic evolution equation

\begin{equation}
\dd \PDMP_\ti = \PDMPDyn(\PDMP_\ti) \dd \ti + [\PDMPTarget(\PDMP_\ti) - \PDMP_\ti] \dd\out_\ti
\end{equation}

and jumps of $\out_\ti$ occur with intensity $\intensityFct(\PDMP_\ti)$. Then for a differentiable initial condition $\PDMPDens(0, \PDMPState)$, the probability density evolves according to the PDE

\begin{align}
&\partial_\ti \PDMPDens(\ti, \PDMPState) = \nonumber \\
&-\partial_\PDMPState(\PDMPDyn(\PDMPState)\PDMPDens(\ti, \PDMPState)) - \intensityFct(\PDMPState)\PDMPDens(\ti, \PDMPState) + \PDMPTargInv'(\PDMPState) \intensityFct(\PDMPTargInv(\PDMPState)) \PDMPDens(\ti, \PDMPTargInv(\PDMPState)).
\label{eq:Poisson-Liouville}
\end{align}

Suppose that the function $\PDMPHelp(\PDMPState)$ evolves as

\begin{equation}
\PDMPHelp'(\PDMPState) = - \frac{\intensityFct(\PDMPState)}{\PDMPDyn(\PDMPState)}\PDMPHelp(\PDMPState)
\end{equation}

with arbitrary initial condition. If $\PDMPDens(\PDMPState)$ fulfills

\begin{equation}
\PDMPDens(\PDMPState) = -\int_{\PDMPTargInv(\PDMPState)}^\UpperBound \frac{\intensityFct(\PDMPStatefrom) \PDMPHelp(\PDMPState)}{\PDMPDyn(\PDMPState)\PDMPHelp(\PDMPTarget(\PDMPStatefrom))} \PDMPDens(\PDMPStatefrom) \dd \PDMPStatefrom,
\label{eq:DirectMethod}
\end{equation}

then the stationarity condition is satisfied:

\begin{equation}
0 = -\partial_\PDMPState(\PDMPDyn(\PDMPState)\PDMPDens( \PDMPState)) - \intensityFct(\PDMPState)\PDMPDens(\PDMPState) + \PDMPTargInv'(\PDMPState) \intensityFct(\PDMPTargInv(\PDMPState)) \PDMPDens( \PDMPTargInv(\PDMPState)).
\label{eq:station}
\end{equation}

Eq. (<ref>) is derived in the Appendix <ref>. Eq. (<ref>) follows from (<ref>) by Leibniz differentiation under the integral sign.

The linear fixed point equation (<ref>) can be used to return a numerical approximation of $\acid(\aci)$ directly. However, it can have a singularity at the equilibrium $\PDMPState$ characterized by $\PDMPDyn(\PDMPState) = 0$. And the ad-hoc discretization is not guaranteed to be a stochastic matrix as in eq. (<ref>).

For the Hawkes process <ref> the equation (<ref>) yields

\begin{equation}
\aciKernel(\aci, \acifrom) = \frac{\acifrom(\aci - \HawkesBase)^{\frac{ \HawkesBase}{\HawkesDecay}-1}}{\HawkesDecay(\acifrom + \HawkesJump - \HawkesBase)^{-\frac{ \HawkesBase}{\HawkesDecay}}} e^{\frac{1}{\HawkesDecay} (\aci -\acifrom - \HawkesJump)}
\label{eq:direct_Hawkes}
\end{equation}

in  (<ref>). By "iterative solution" <cit.> this fixed point iteration could have originally been meant instead of the method of steps.
For <ref> note that eq. (<ref>) corrects eq.<cit.> by a missing factor. The equation (<ref>) yields

\begin{equation}
\aciKernel(\aci, \acifrom) = \frac{\acifrom(\aci - \RTRoot_1)^{\frac{\RTRoot_1}{\Delta \RTRoot} - 1}(\RTRoot_2 - \PDMPTargInv(\acifrom))^\frac{\RTRoot_2}{\Delta \RTRoot}}{(\PDMPTargInv(\acifrom) - \RTRoot_1)^{\frac{\RTRoot_1}{\Delta \RTRoot}}(\RTRoot_2 - \aci)^{\frac{\RTRoot_2}{\Delta \RTRoot} + 1}}.
\label{eq:direct_leakage}
\end{equation}

Fig. <ref> (Appendix <ref>) shows agreement for three Hawkes examples, while in fig. <ref> the accuracy for the considered random telegraph with dark current suffers, possibly caused by the singularity at the equilibrium $\RTRoot_1$.

scale only axis,
xlabel = conditional intensity $\hat{\lambda}_\infty$,
ylabel = ACID,
axis background/.style=fill=white,
legend style = at =(0,1), anchor = north west, legend cell align=left,align = left
[ybar interval, fill=mycolor0, fill opacity=0.3, draw=gray, area legend] table[row sep=crcr]
x	y

0	9.09090909090909e-05

0.011	0

0.022	0

0.033	0

0.044	0

0.055	0

0.066	0

0.077	0

0.088	0

0.099	0

0.11	0

0.121	0

0.132	0

0.143	0

0.154	0

0.165	0

0.176	0

0.187	18.4201818181818

0.198	3.31818181818182

0.209	1.99590909090909

0.22	1.51818181818182

0.231	1.25827272727273

0.242	1.054

0.253	0.94690909090909

0.264	0.864363636363636

0.275	0.805363636363636

0.286	0.744181818181818

0.297	0.711454545454545

0.308	0.671363636363639

0.319	0.634999999999999

0.33	0.60990909090909

0.341	0.59090909090909

0.352	0.575727272727272

0.363	0.571454545454545

0.374	0.56390909090909

0.385	0.545545454545457

0.396	0.544636363636363

0.407	0.541090909090909

0.418	0.522818181818181

0.429	0.527272727272727

0.44	0.532727272727272

0.451	0.523818181818184

0.462	0.52290909090909

0.473	0.528272727272727

0.484	0.515454545454545

0.495	0.52818181818182

0.506	0.520818181818181

0.517	0.530999999999999

0.528	0.512636363636363

0.539	0.546545454545454

0.55	0.541818181818181

0.561	0.564454545454545

0.572	0.562272727272727

0.583	0.575454545454545

0.594	0.591181818181818

0.605	0.586454545454545

0.616	0.55390909090909

0.627	0.481000000000004

0.638	0.475454545454545

0.649	0.467272727272727

0.66	0.475909090909091

0.671	0.484545454545454

0.682	0.485818181818181

0.693	0.498090909090909

0.704	0.495

0.715	0.509545454545454

0.726	0.510181818181818

0.737	0.53090909090909

0.748	0.559636363636363

0.759	0.560363636363636

0.77	0.584000000000005

0.781	0.584818181818181

0.792	0.623818181818181

0.803	0.632727272727272

0.814	0.647181818181818

0.825	0.683181818181818

0.836	0.695727272727272

0.847	0.739636363636363

0.858	0.756454545454545

0.869	0.803818181818181

0.88	0.844181818181817

0.891	0.880545454545454

0.902	0.928999999999999

0.913	0.979090909090918

0.924	1.03

0.935	1.08981818181818

0.946	1.17736363636364

0.957	1.247

0.968	1.33145454545454

0.979	1.42681818181818

0.99	1.54109090909091

1.001	1.62072727272729

1.012	1.72481818181816

1.023	1.86063636363638

1.034	2.02981818181816

1.045	2.19018181818184

1.056	2.3452727272727

1.067	2.39245454545457

1.078	2.15763636363634

1.089	1.05145454545456

1.1	1.05145454545456

Monte Carlo samples
[color=mycolor1, line width = 1.5pt]
table[row sep=crcr]
0.190325524128568	83.9275919518654

0.190780475104248	34.3673300862448

0.191235426079928	22.7025725820388

0.192145328031273	14.1005338040288

0.193055229982633	10.4749818709088

0.193965131933993	8.43752672500014

0.194875033885353	7.11872472797033

0.196239886812393	5.81943238931468

0.197604739739432	4.95741401911877

0.198969592666472	4.34019803579127

0.200789396569192	3.74362428326594

0.202609200471912	3.30724490981362

0.204429004374632	2.9730305843086

0.206248808277351	2.70819247662064

0.208523563155751	2.4449773480496

0.210798318034151	2.23512575363308

0.21307307291255	2.0636063044077

0.215347827790936	1.92060269309785

0.218077533645015	1.77736061153891

0.220807239499095	1.65767734730275

0.223536945353175	1.55608956307942

0.226266651207254	1.46871977801769

0.229451308037014	1.38102277600595

0.232635964866773	1.30555245666957

0.235820621696519	1.23989217661943

0.239005278526278	1.18223009802476

0.242644886331718	1.12436078520594

0.246284494137157	1.07354941846332

0.249924101942597	1.02857912205448

0.253563709748036	0.988501153287075

0.257658268529141	0.948332178586114

0.261752827310261	0.912560335368354

0.26584738609138	0.880513554700926

0.2699419448725	0.851652210713226

0.274491454629299	0.822789513496303

0.279040964386084	0.796830625365047

0.283590474142883	0.77337772969075

0.288139983899683	0.752104004013191

0.292689493656482	0.732738377924989

0.297693954388961	0.713370484919579

0.302698415121426	0.695785112078354

0.307702875853906	0.679771051945238

0.312707336586385	0.665150335980812

0.317711797318864	0.651771927784111

0.323171209027009	0.638443041973986

0.328630620735169	0.626296329802429

0.334090032443328	0.61520973590828

0.339549444151487	0.605078510168113

0.345008855859632	0.595812261515817

0.350468267567791	0.587332598035431

0.35638263025163	0.578954933994126

0.36229699293547	0.571344810518113

0.368211355619295	0.564436528704121

0.374125718303134	0.558172754848457

0.380040080986973	0.552503269988776

0.386409394646492	0.547011786479089

0.392778708305997	0.542110062988399

0.399148021965516	0.53775574978188

0.405517335625035	0.533911673005818

0.41188664928454	0.530545132361553

0.418710913919739	0.527435378493067

0.425535178554938	0.524809652912168

0.432359443190123	0.522640712448265

0.439183707825322	0.520904717168818

0.446462923436201	0.519506752469098

0.453742139047066	0.518555703257519

0.461021354657944	0.518032598896312

0.468300570268823	0.517921154573159

0.476034736855368	0.518238311226341

0.483768903441927	0.51899055813729

0.491503070028486	0.520166389440504

0.499237236615031	0.521756553825313

0.507426354177269	0.523883911389291

0.515615471739508	0.526461542625725

0.523804589301733	0.529485472623648

0.531993706863972	0.532953817795686

0.540637775401891	0.537097124317782

0.549281843939795	0.541738143035531

0.557925912477714	0.546881573245997

0.566569981015618	0.552534351514126

0.575214049553537	0.558705630280755

0.583858118091456	0.565406781088996

0.59250218662936	0.572651421292548

0.601146255167279	0.580455463618875

0.609790323705184	0.588837188421863

0.618434392243103	0.597817338892739

0.621619049072862	0.601281222902699

0.622074000048542	0.572512940983188

0.622983901999902	0.573472728922681

0.623438852975582	0.561951098038392

0.624348754926942	0.562900651912244

0.624803705902622	0.555435901951697

0.625713607853982	0.556381834276849

0.626168558829661	0.55080175871224

0.627078460781021	0.551747126811065

0.627533411756701	0.547273790359313

0.628443313708047	0.548220386620159

0.628898264683727	0.544482100596696

0.629808166635087	0.545431117120174

0.630263117610767	0.542219877693043

0.631173019562127	0.543172173034591

0.631627970537807	0.540359186985825

0.632537872489166	0.541315414732864

0.632992823464846	0.53881509770693

0.633902725416206	0.539775776344271

0.634357676391886	0.53752820583685

0.635267578343246	0.538493760080229

0.635722529318926	0.536455256821682

0.636632431270286	0.537426043681023

0.637087382245966	0.535563732943885

0.637997284197326	0.53654005899466

0.638452235173006	0.534828544631296

0.639362137124365	0.53581067792507

0.639817088100045	0.534229911978358

0.640726990051405	0.535218090529582

0.641181941027085	0.533751955469029

0.641636892002765	0.534248348119263

0.642091842978445	0.532883766325597

0.643001744929805	0.533881344536525

0.643456695905485	0.532607262613553

0.644366597856845	0.533611461020058

0.644821548832525	0.532418763113029

0.645731450783884	0.53342974756977

0.646186401759564	0.532310719739826

0.647096303710924	0.533328645502039

0.647551254686604	0.532276699840821

0.648006205662284	0.532788364298725

0.648461156637964	0.531797591943075

0.649371058589324	0.532826450905972

0.649826009565004	0.531892118001139

0.65073591151635	0.532928319359229

0.65119086249203	0.532046132429841

0.652100764443389	0.533089802974004

0.652555715419069	0.532255990212533

0.653010666394749	0.532780763876886

0.653465617370429	0.531991577419049

0.654375519321789	0.533047105893118

0.654830470297469	0.532299965466024

0.655740372248829	0.533363299997418

0.656195323224509	0.532655569615073

0.656650274200189	0.533190326448462

0.657105225175869	0.532519119627921

0.658015127127229	0.533594906579893

0.658470078102908	0.532958592953463

0.659379980054268	0.534042494587027

0.659834931029948	0.533439160103896

0.660289882005628	0.533984338463455

0.660744832981308	0.533411629632553

0.661654734932668	0.534508528848718

0.662109685908348	0.533965462924897

0.663019587859708	0.535070767522072

0.663474538835388	0.534555912163725

0.663929489811068	0.535111921320848

0.664384440786748	0.534623263412513

0.665294342738108	0.535742078433543

0.665749293713787	0.535279098763709

0.666659195665147	0.536406598815503

0.667114146640827	0.535968199347309

0.667569097616507	0.536535429958732

0.668024048592187	0.536119818963158

0.668933950543547	0.537261321539944

0.669388901519227	0.536868325728648

0.669843852494907	0.537442634828764

0.670298803470587	0.537070566457189

0.671208705421947	0.538226379134684

0.671663656397627	0.537875249900978

0.672118607373292	0.538456790004716

0.672573558348972	0.538124990454961

0.673483460300332	0.539295416143759

0.673938411276012	0.538983116061743

0.674393362251692	0.539572037757509

0.674848313227372	0.539277686117245

0.675758215178732	0.540463024615136

0.676213166154412	0.540186928140997

0.676668117130092	0.540783380799311

0.677123068105772	0.540524031331557

0.678032970057131	0.541724580621747

0.678487921032811	0.541482405022734

0.678942872008491	0.542086537392663

0.679397822984171	0.541860058023573

0.680307724935531	0.543076115374177

0.680762675911211	0.542865865755772

0.681217626886891	0.543477826473193

0.681672577862571	0.543282349577382

0.682582479813931	0.544514212484103

0.683037430789611	0.54433413758656

0.683492381765291	0.544954075612551

0.683947332740971	0.544787958250055

0.684402283716651	0.545410949965813

0.684857234692331	0.545258409421862

0.68576713664369	0.546512546578114

0.68622208761937	0.546374338079204

0.68667703859505	0.54700551957967

0.68713198957073	0.546880207204637

0.68804189152209	0.548150874603024

0.68849684249777	0.54803928992601

0.68895179347345	0.548678812186594

0.68940674444913	0.548579517443684

0.68986169542481	0.549222246001918

0.69031664640049	0.549134954558411

0.69122654835185	0.55042893159586

0.69168149932753	0.550354577845354

0.692136450303209	0.551005863172406

0.692591401278889	0.55094300076216

0.693046352254569	0.551597582756969

0.693501303230249	0.551545967821269

0.694411205181595	0.552863870179621

0.694866156157275	0.552824530372746

0.695321107132955	0.553487888196926

0.695776058108635	0.55345936126966

0.696231009084315	0.554126107073074

0.696685960059995	0.554108184450598

0.697595862011354	0.555450636532015

0.698050812987034	0.555444431844649

0.698505763962714	0.556120176308497

0.698960714938394	0.556124205552905

0.699415665914074	0.556803430571847

0.699870616889754	0.556817514597682

0.700780518841114	0.558185151523816

0.701235469816794	0.558210484380353

0.701690420792474	0.558898935299666

0.702145371768154	0.558934006424025

0.702600322743834	0.559626032024781

0.703055273719514	0.559670687036657

0.703510224695194	0.560366315136363

0.703965175670874	0.560420407053314

0.704875077622233	0.561821149308415

0.705330028597913	0.561885971429405

0.705784979573593	0.562591125773864

0.706239930549273	0.562665127529513

0.706694881524953	0.563373980803121

0.707149832500633	0.563457035167502

0.707604783476313	0.564169615817178

0.708059734451993	0.564261601763192

0.708514685427673	0.564977938376771

0.708969636403353	0.56507874058093

0.709879538354713	0.566521278185746

0.710334489330393	0.5666323053631

0.710789440306073	0.56735854816516

0.711244391281753	0.567478193107618

0.711699342257432	0.568208291616628

0.712154293233112	0.56833645454526

0.712609244208792	0.569070438125905

0.713064195184472	0.569207023966527

0.713519146160152	0.569944922157745

0.713974097135832	0.570089840327697

0.714429048111512	0.570831682848706

0.714883999087192	0.570984847056565

0.715793901038538	0.572478863389918

0.716248852014218	0.572641812958651

0.716703802989898	0.573394030752127

0.717158753965577	0.573565082319121

0.717613704941257	0.574321348776849

0.718068655916937	0.574500427889319

0.718523606892617	0.575260773663629

0.718978557868297	0.575447809475833

0.719433508843977	0.576212265428907

0.719888459819657	0.576407190554207

0.720343410795337	0.577175787763466

0.720798361771017	0.577378538134965

0.721253312746697	0.578151307898423

0.721708263722377	0.578361822637433

0.722163214698057	0.579138796479

0.722618165673737	0.57935701777096

0.723073116649417	0.580138227445588

0.723528067625097	0.580364100422813

0.723983018600777	0.581149577921721

0.724437969576456	0.58138305055266

0.724892920552136	0.582172828108412

0.725347871527816	0.582413851092937

0.725802822503496	0.583207961184485

0.726257773479176	0.583456487854875

0.726712724454856	0.584254963212658

0.727167675430536	0.584510949439832

0.728077577381896	0.586119281457371

0.728532528357576	0.586384532267886

0.728987479333256	0.587194438239194

0.729442430308936	0.58746708506564

0.729897381284616	0.5882814722886

0.730352332260296	0.588561478203147

0.730807283235976	0.58938038063765

0.731262234211655	0.589667710926634

0.731717185187335	0.590491162810466

0.732172136163015	0.590785784902806

0.732627087138695	0.591613820756521

0.733082038114375	0.591915704155838

0.733536989090055	0.59274835878773

0.734446891041415	0.593371826454671

0.734901842017095	0.594211106023096

0.735356792992775	0.594532620895947

0.735811743968455	0.595376607953952

0.736266694944135	0.595705263321776

0.736721645919815	0.596553993691643

0.737176596895495	0.596889768403841

0.73763154787116	0.597743278218374

0.73808649884684	0.598086152855217

0.73854144982252	0.598944478562615

0.7389964007982	0.599294435387932

0.73945135177388	0.600157613756792

0.73990630274956	0.600514636673182

0.74036125372524	0.601382704797402

0.74081620470092	0.601746779303426

0.7412711556766	0.602619774607234

0.74172610665228	0.602990887756945

0.74218105762796	0.603868847999863

0.74263600860364	0.604246988364167

0.74309095957932	0.60512995164602

0.743545910555	0.605515109276084

0.744000861530679	0.606403114042095

0.744455812506359	0.606795280434383

0.744910763482039	0.607688365480215

0.745820665433399	0.608491488130127

0.746275616409079	0.60939190204364

0.746730567384759	0.609802828309981

0.747185518360439	0.610708421087068

0.747640469336119	0.611126315835804

0.748095420311799	0.612037127513673

0.748550371287479	0.612461988839726

0.749005322263159	0.613378059829671

0.749460273238839	0.613809887094291

0.750825126165878	0.61609676436413

0.751735028117238	0.616992917156281

0.753099881044278	0.619324584445749

0.753554832019958	0.619788889414053

0.754919684946998	0.622156436265698

0.755829586898358	0.623117467514277

0.757194439825398	0.625530700245505

0.757649390801078	0.626027437065915

0.759014243728117	0.628477380872809

0.759924145679463	0.629503324531782

0.761288998606503	0.631999962988147

0.761743949582183	0.632529208796385

0.762653851533543	0.634061516339187

0.763108802509223	0.634602241427302

0.764473655436262	0.6371651832563

0.764928606411942	0.637720035227773

0.765838508363302	0.639298620280201

0.766293459338982	0.639864974196044

0.767658312266022	0.642495513875943

0.768568214217382	0.643661104522224

0.769933067144422	0.646340791054044

0.770842969095781	0.647543831246466

0.772207822022821	0.650273383214937

0.773117723974181	0.651514068609714

0.774482576901221	0.654294238816064

0.775392478852581	0.655572801685409

0.776757331779621	0.658404377117776

0.77766723373098	0.659721085252116

0.77857713568234	0.661493543382377

0.77903208665802	0.662165862772454

0.78039693958506	0.665081371731333

0.781306841536406	0.666460259546156

0.782216743487766	0.668290866592571

0.782671694463446	0.668994419683912

0.784036547390485	0.671996065402709

0.784946449341845	0.673437971104477

0.785856351293205	0.675328089678175

0.786311302268885	0.676063329008329

0.787221204220245	0.677976182834215

0.787676155195925	0.678723571342175

0.789041008122965	0.681847468892599

0.789950910074324	0.683377913658759

0.790860812025684	0.685352319282032

0.791315763001364	0.686132119425892

0.792225664952724	0.688130049165366

0.793135566904084	0.690143714188935

0.794045468855444	0.691748706436769

0.794955370806804	0.693794368616864

0.796320223733844	0.696681598060749

0.797685076660883	0.699605983789198

0.799504880563603	0.703423577375318

0.801324684466323	0.707309611924643

0.802234586417683	0.709487358705033

0.804509341296068	0.714381099250332

0.805419243247428	0.716618963800784

0.807239047150148	0.720731963393419

0.808148949101508	0.72302259577485

0.810878654955587	0.729174988080061

0.811788556906947	0.731537731393558

0.814063311785347	0.73690763022482

0.814973213736707	0.739335309087949

0.817702919590786	0.745911284937705

0.818612821542146	0.748415275292942

0.821342527396226	0.755226359513998

0.822252429347586	0.757809093634151

0.825437086177331	0.76598770184718

0.826346988128691	0.768661997834513

0.82998659593413	0.778347828461051

0.83089649788549	0.781127790306428

0.83453610569093	0.791247155668614

0.837720762520689	0.800446751430712

0.839085615447729	0.80470735087188

0.843180174228849	0.816995432197999

0.847274733009954	0.829774147731229

0.851824242766753	0.844525638797833

0.853189095693793	0.849362928911418

0.858193556426272	0.866122883096779

0.859558409353312	0.871247262172716

0.866382673988511	0.895793040147154

0.873206938623696	0.922035380856201

0.881396056185935	0.955519940838727

0.891404977650879	0.999579283413667

0.901868850091518	1.05032348168938

0.915517379361901	1.12362146453698

0.920521840094381	1.15295063824713

0.9223416439971	1.1637574246695

0.92552630082686	1.18343418970298

0.92734610472958	1.19485225772364

0.930530761559339	1.21563869075666

0.933260467413419	1.23334122488171

0.935990173267484	1.25209525865854

0.938719879121564	1.27092495427559

0.941449584975643	1.29086116571183

0.944634241805403	1.31413041001944

0.947818898635163	1.33826629857734

0.951458506440602	1.36680576542926

0.955098114246042	1.39658302065084

0.960102574978507	1.43866369644761

0.965107035710986	1.48332589548386

0.971931300346185	1.547266818163

0.980575368884089	1.63447811833888

0.992404094251768	1.76596683445199

0.995588751081527	1.80402209606525

0.998773407911287	1.84282680583642

1.00150311376535	1.87757358982348

1.00650757449783	1.9421243054061

1.01242193718167	2.0226064128743

1.02425066254933	2.19266087933553

1.02789027035477	2.250893293236

1.03107492718453	2.30428068845438

1.03653433889269	2.39805642899664

1.05609723084692	2.74465053122614

1.05700713279828	2.75907117286052

1.05928188767668	2.79567369354984

1.06155664255508	2.82766840102209

1.06246654450644	2.8412880216046

1.06610615231186	2.88434314499678

1.06701605426322	2.89404415540419

1.0674710052389	2.89630207785751

1.06838090719026	2.90514420566844

1.06883585816594	2.9068206526118

1.0697457601173	2.91468095387587

1.07020071109298	2.91570179449573

1.07065566206866	2.91916342398164

1.07156556402002	2.91991276320942

1.0720205149957	2.92270290874906

1.07247546597138	2.9223640495233

1.07338536792274	2.92069152312294

1.0742952698741	2.91758892554559

1.07520517182546	2.91293390573361

1.0765700247525	2.9027495072278

1.07702497572818	2.89842750659707

1.07793487767954	2.88453267080145

1.07838982865522	2.87857521847138

1.0811195345093	2.81248091478724

1.08293933841202	2.74467694870233

1.08475914231474	2.65856422528553

1.0874888481688	2.46721834988581

1.0897636030472	2.24111801166954

1.09249330890128	1.84766127012635

1.09476806377968	1.38375328435495

1.0965878676824	0.907001905389762

1.09931757353648	0.0745892278449531

1.09977252451216	0.00557725417623089

;direct method
[color=mycolor3, line width = 1.5pt, dashed]
table[row sep=crcr]
0.199879494617834	4.05133135193924

0.201244347544873	3.6535529967389

0.203064151447591	3.24326798913469

0.20488395535031	2.92608549330264

0.206703759253029	2.67294448897628