Coarse-Grained Hawkes Processes

1. Introduction

Generally, past events in natural and social systems facilitate the occurrence of future events through a self-exciting mechanism. The Hawkes process [1,2], a class of self-exciting point processes, has become a popular tool for modeling such processes continuously. Its defining intensity function, or instantaneous event rate, is composed of a baseline rate augmented by the cumulative influence of prior events, thereby capturing the self-exciting mechanism wherein each event elevates the probability of subsequent occurrences. Owing to the prevalence of self-excitation across diverse domains, the Hawkes process has been extensively employed in a wide range of disciplines, including seismology [3,4], neurophysiology [5,6], genomics [7], finance [8,9], social media analytics [10,11,12], criminology [13,14], terrorism studies [15], and traffic incident analysis [16].

The Hawkes process is applicable to scenarios wherein individual events are distinguishable, as it characterizes a series of discrete occurrences along a continuous timeline. Nevertheless, it proves insufficient for modeling count-based sequences in which the precise timing of events is unobserved and only their aggregated counts within successive intervals are available. A representative case arises in epidemiological studies, where individual infection events are not monitored in real time, but rather, the daily incidence rates are recorded. Consequently, there is a need to develop a count-based time series model that retains the salient features of the Hawkes process for such contexts.

Multiple methodologies have been proposed for fitting the Hawkes process to bin-count data. Kirchner demonstrated that the distribution of bin-count sequences generated by Hawkes processes can be effectively approximated using an integer-valued autoregressive (INAR) model, wherein conditional least squares estimation is employed to infer the underlying Hawkes process [17,18]. Shlomovich et al. introduced a Monte Carlo expectation-maximization (MC-EM) framework, which integrates an efficient sampling algorithm for the latent event times with an EM procedure to maximize the likelihood function [19,20]. Similarly, Chen et al. developed a Pseudo-Marginal Metropolis-Hastings (PMMH) algorithm for maximum likelihood estimation, wherein the likelihood is approximated via a sequential Monte Carlo approach [21]. Alternatively, Cheysson and Lang advocated for a spectral estimation technique grounded in the Whittle likelihood within the univariate context [22].

Although Kirchner’s approach provides consistent and asymptotically normal estimators for the underlying Hawkes process as the bin width approaches zero, it introduces bias when the bin size exceeds the typical support of the excitation function, due to the INAR model’s neglect of intra-bin excitation dynamics. In contrast, the MC-EM, PMMH and spectral estimation methods yield less biased results by accounting for excitation effects occurring within each bin. Nonetheless, the MC-EM and PMMH techniques are computationally intensive, whereas the spectral method is validated only in the univariate setting.

Accordingly, the objective of this study is to develop a computationally efficient methodology that extends naturally to the multivariate setting. Our proposed approach involves constructing a count-based time series model that approximates the bin-count sequence generated by Hawkes processes, upon which an estimation procedure is built. Unlike Kirchner’s method, which relies on a straightforward discretization, namely, evaluating the intensity function at discrete time points, to relate the Hawkes and INAR processes [17], our model is derived through a coarse graining procedure. The resulting framework, termed the coarse-grained Hawkes process, is defined in a conceptually simple and analytically tractable form, while effectively preserving the second-order statistical characteristics of the bin-count Hawkes process, even in regimes where the bin size is large relative to the excitation function’s effective range.

The structure of this paper is as follows. Section 2 provides a concise overview of Hawkes processes. Section 3 introduces the coarse-grained Hawkes process along with the proposed estimation methodology. In Section 4, we assess the approximation accuracy of the coarse-grained Hawkes process in representing the bin-count Hawkes process. Additionally, a simulation study is conducted to benchmark the performance of the proposed estimation technique against existing approaches. Section 5 concludes with a discussion of the findings.

2. Review of Hawkes Processes

In this section, we present a brief review of Hawkes processes. Let ${\left\{t_i\right\}}_{i\in \mathbb{N}}$ denote a sequence of non-negative random variables representing the occurrence times of events on ${\mathbb{R}}_{+}$, satisfying $t_i<t_{i+1}$ for $i\in \mathbb{N}$. Define the associated counting process by $N\left(t\right)={\sum }_{i\in \mathbb{N}}{\mathbf{1}}_{t_i\le t}$, and let ${\mathcal{H}}_t=\{N\left(u\right):u<t\}$ denote the history of events up to, but not including, time t. We consider a point process such that

$$\begin{array}{cc}\hfill P(N(t+\delta t)-N\left(t\right)=1\mid {\mathcal{H}}_t)& =\lambda \left(t\right)\delta t+o\left(\delta t\right),\hfill \\ \hfill P(N(t+\delta t)-N\left(t\right)>1\mid {\mathcal{H}}_t)& =o\left(\delta t\right),\hfill \end{array}$$

as $\delta t\to 0$, where $\lambda \left(t\right)$ is the conditional intensity function that uniquely characterizes the point process. A univariate Hawkes process is then defined by the conditional intensity function

$$\lambda \left(t\right)=\mu +{\int }_0^t\varphi (t-u)dN\left(u\right),$$

where $\mu \ge 0$ is the baseline intensity and $\varphi (·)$ is a non-negative excitation kernel satisfying $\varphi \left(t\right)=0$ for $t<0$ [1]. This formulation captures the self-exciting nature of the process, wherein the intensity at any given time depends on the historical sequence of events. The branching ratio is defined as

$$\alpha :={\int }_0^{\infty }\varphi \left(t\right)dt,$$

which quantifies the expected number of subsequent events triggered by a single occurrence. The process admits a stationary distribution provided that $\alpha <1$.

The univariate Hawkes process naturally generalizes to the multivariate case [1,2]. For $d\in \mathbb{N}$, a d-dimensional Hawkes process $\mathbf{N}\left(t\right)={(N_1\left(t\right),\dots ,N_d\left(t\right))}^{\mathsf{T}}$ comprises d jointly defined point processes on ${\mathbb{R}}_{+}$. The corresponding vector-valued conditional intensity function is given by

$$\mathbf{\lambda }\left(t\right)=\mathbf{\mu }+{\int }_0^t\Phi (t-s)d\mathbf{N}\left(s\right),$$

where $\mathbf{\mu }={({\mu }_1,\dots ,{\mu }_d)}^{\mathsf{T}}\in {\mathbb{R}}_{\ge 0}^d$ is the baseline intensity vector, and the excitation kernel $\Phi ={\left({\varphi }_{ij}\right)}_{1\le i,j\le d}$ is a matrix-valued function with non-negative entries satisfying ${\varphi }_{ij}\left(t\right)=0$ for $t<0$. Each entry ${\varphi }_{ij}(·)$ describes the influence of process j on the intensity of process i, thus capturing both self-excitation and mutual excitation among the components. The multivariate Hawkes process is asymptotically stationary if the spectral radius of the branching matrix

$$A:={\int }_0^{\infty }\Phi \left(t\right)dt$$

is strictly less than one.

Finally, we briefly recall the martingale properties of point processes that are instrumental for our methodological developments. For a comprehensive treatment, see [23]. Define the process $\mathbf{M}\left(t\right)=\mathbf{N}\left(t\right)-\Lambda \left(t\right)$, where $\Lambda \left(t\right)={\int }_0^t\mathbf{\lambda }\left(u\right)du$. Then $\mathbf{M}\left(t\right)$ is a martingale with respect to the filtration ${\mathcal{H}}_t$, satisfying the property that the conditional expectation of its increment is zero:

$$E\left[\mathbf{M}\left(t\right)-\mathbf{M}\left(s\right)\mid {\mathcal{H}}_s\right]=\mathbf{0},t>s.$$

Moreover, the conditional covariance matrix of the martingale increment is given by

$$E\left[(\mathbf{M}\left(t\right)-\mathbf{M}\left(s\right)){(\mathbf{M}\left(t\right)-\mathbf{M}\left(s\right))}^{\mathsf{T}}\mid {\mathcal{H}}_s\right]=\mathrm{diag}\left(E[\Lambda \left(t\right)-\Lambda \left(s\right)\mid {\mathcal{H}}_s]\right),t>s,$$

which follows from the quadratic variation of the martingale.

3. Coarse-Grained Hawkes Process

3.1. Motivation

We consider a scenario wherein the exact timing of individual events is unobserved; instead, we observe an aggregated count of these latent continuous-time events within discrete time intervals. We define

$${\mathbf{X}}_n=\mathbf{N}\left(n\Delta t\right)-\mathbf{N}\left((n-1)\Delta t\right),n\in \mathbb{N},$$

as the bin-wise event counts over intervals of size $\Delta t$. When $\{\mathbf{N}\left(t\right):t\in {\mathbb{R}}_{+}\}$ constitutes a Hawkes process, the corresponding discrete-time vector process $\{{\mathbf{X}}_n:n\in \mathbb{N}\}$ is referred to as the binned Hawkes process [19,20].

Although a closed-form representation of the probability distribution of binned Hawkes processes is not available, an approximate formulation can be derived by discretizing the conditional intensity function as follows:

$${\mathbf{\lambda }}_n={({\lambda }_{1n},\dots ,{\lambda }_{dn})}^{\mathsf{T}}=\mathbf{\mu }\Delta t+\sum _{k=1}^{n-1}\Phi \left(k\Delta t\right)\Delta t{\mathbf{X}}_k,$$

(1)

and assuming that the bin counts follow Poisson distributions:

$${\mathbf{X}}_n\mid {\mathbf{X}}_1,\dots ,{\mathbf{X}}_{n-1}\sim \prod _{j=1}^d\mathrm{Poisson}\left({\lambda }_{jn}\right).$$

(2)

This method, known as the binned Poisson approximation, converges in distribution to the true Hawkes process as $\Delta \to 0$. However, its accuracy degrades with larger $\Delta t$ due to discretization errors in (1) and the conditional independence assumption in (2), which neglects the intra-bin excitation effect. Nonetheless, it facilitates a tractable estimation procedure: the log-likelihood of the sequence ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_n$ is given by

$$logP({\mathbf{X}}_1,\dots ,{\mathbf{X}}_n)=\sum _{k=1}^n\sum _{j=1}^d\left(X_{jk}log{\lambda }_{jk}-{\lambda }_{jk}\right).$$

Parameter estimation is then performed via maximization of this log-likelihood function [19,20,24].

In light of the limitations of this approximation for larger $\Delta t$, we propose an alternative count time series model that more accurately approximates the binned Hawkes process while retaining computational simplicity for inference. We elucidate the underlying heuristic in the univariate setting.

The central idea is to replace the crude discretization with the expected value of the conditional intensity integrated over each time bin. The expected count within the nth bin is ${\int }_{(n-1)\Delta t}^{n\Delta t}\lambda \left(t\right)dt$, which depends on the trajectory $\left\{N\right(u):0<u\le n\Delta t\}$. Consequently, we consider the conditional expectation, given the bin counts $\{X_1,\dots ,X_n\}$:

$$\begin{array}{cc}\hfill {\lambda }_n& =E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\lambda \left(t\right)dt\mid X_1,\dots ,X_n\right]\hfill \\ \hfill & =\sum _{k=1}^{n-1}E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\left({\int }_{(k-1)\Delta t}^{k\Delta t}\varphi (t-u)dN\left(u\right)\right)dt\mid X_k\right]\hfill \\ \hfill & +E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\left({\int }_{(n-1)\Delta t}^t\varphi (t-u)dN\left(u\right)\right)dt\mid X_n\right].\hfill \end{array}$$

(3)

As these conditional expectations do not admit closed-form expressions, we approximate them under the assumption that event times are uniformly distributed within each bin. Accordingly, the first term on the right-hand side of (3) is approximated as

$$\begin{array}{cc}\hfill & E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\left({\int }_{(k-1)\Delta t}^{k\Delta t}\varphi (t-u)dN\left(u\right)\right)dt\mid X_k\right]\hfill \\ \hfill & \approx {\int }_{(k-1)\Delta t}^{k\Delta t}\cdots {\int }_{(k-1)\Delta t}^{k\Delta t}{\displaystyle \frac{1}{\Delta t^{X_k}}}\left(\sum _{i=1}^{X_k}{\int }_{(n-1)\Delta t}^{n\Delta t}\varphi (t-t_i)dt\right)d{t}_1\cdots d{t}_{X_k}\hfill \\ \hfill & ={\varphi }_{n-k}{X}_k,\hfill \end{array}$$

where ${\varphi }_{n-k}$ denotes the coarse-grained kernel, defined as

$${\varphi }_{n-k}={\displaystyle \frac{1}{\Delta t}}{\int }_{(k-1)\Delta t}^{k\Delta t}{\int }_{(n-1)\Delta t}^{n\Delta t}\varphi (t-u)dudt.$$

(4)

Similarly, the second term on the right-hand side of (3) is approximated by

$$\begin{array}{cc}\hfill & E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\left({\int }_{(n-1)\Delta t}^t\varphi (t-u)dN\left(u\right)\right)dt\mid X_n\right]\hfill \\ \hfill & =E\left[{\int }_{(n-1)\Delta t}^{n\Delta t}\left({\int }_{(n-1)\Delta t}^{n\Delta t}\varphi (t-u)dN\left(u\right)\right)dt\mid X_n\right]\approx {\varphi }_0{X}_n,\hfill \end{array}$$

where we have utilized the causal property $\varphi \left(t\right)=0$ for $t<0$. Substituting these approximations back into (3) yields

$${\lambda }_n\approx \mu \Delta t+\sum _{k=1}^n{\varphi }_{n-k}{X}_k.$$

(5)

In contrast to the formulation in (1), equation (5) incorporates the coarse-grained kernel, capturing both inter-bin excitation and the intra-bin self-excitation effect through the term ${\varphi }_0{X}_n$. Building upon this approximation, we will formally define the coarse-grained Hawkes process in the subsequent section.

3.2. Definition

Firstly, we formally define the coarse-grained kernel, previously derived heuristically in (4).

Definition 1.

Let $\varphi \left(t\right)$ be a nonnegative excitation kernel defined on the real line, satisfying $\varphi \left(t\right)=0$ for $t<0$. The coarse-grained kernel ${\left\{{\varphi }_k\right\}}_{k\in {\mathbb{N}}_0}$ with bin size $\Delta t>0$ is defined by

$${\varphi }_k=\left\{\begin{array}{cc}{\psi }_0,\hfill & k=0,\hfill \\ {\psi }_k-{\psi }_{k-1},\hfill & k=1,2,\dots ,\hfill \end{array}\right.$$

(6)

where

$${\psi }_k={\displaystyle \frac{1}{\Delta t}}{\int }_{k\Delta t}^{(k+1)\Delta t}\left({\int }_0^t\varphi \left(u\right)du\right)dt.$$

(7)

It can be readily verified that the formulation (6)–(7) coincides with the expression in (4).

Lemma 1.

The coarse-grained kernel satisfies

$$\begin{array}{c}\hfill \sum _{k=0}^{\infty }{\varphi }_k={\int }_0^{\infty }\varphi \left(u\right)du=:\alpha .\end{array}$$

Moreover, it holds that ${\varphi }_0<\alpha$.

Proof. 

See Appendix A.1.    □

Lemma 1 guarantees that the total mass of the coarse-grained kernel equals the integral of the excitation kernel, i.e., the branching ratio. Furthermore, it ensures that the kernel does not collapse to the degenerate case ${\varphi }_0=\alpha$ as long as $\Delta t<\infty$.

Based on the coarse-grained kernel, we introduce the coarse-grained Hawkes process. Consider a probability space equipped with a sequence of d-dimensional, integer-valued random vectors ${\left\{{\mathbf{X}}_k\right\}}_{k\in \mathbb{Z}}$. Define the history up to bin k as $H_k=\{{\mathbf{X}}_j:j\le k\}$. We consider the d-variate extension of the coarse-grained process given in (5):

$${\mathbf{\lambda }}_n=\mathbf{\mu }\Delta t+\sum _{k\le n}{\Phi }_{n-k}{\mathbf{X}}_k,n\in \mathbb{Z},$$

(8)

where $\mathbf{\mu }\in {\mathbb{R}}_{\ge 0}^d$ is the baseline intensity vector, and each ${\Phi }_{n-k}\in {\mathbb{R}}^{d\times d}$ is a matrix whose elements are given by the coarse-grained kernels for bin size $\Delta t$. Define the residual process as

$$\Delta {\mathbf{M}}_n={\mathbf{X}}_n-{\mathbf{\lambda }}_n,n\in \mathbb{Z},$$

(9)

and impose the following conditional moment properties, analogous to the martingale conditions in point process theory:

$$E[\Delta {\mathbf{M}}_n\mid H_{n-1}]=\mathbf{0},$$

(10)

$$E[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_n^{\mathsf{T}}\mid H_{n-1}]=\mathrm{diag}\left(E[{\mathbf{\lambda }}_n\mid H_{n-1}]\right).$$

(11)

The coarse-grained Hawkes process is then defined as a sequence $\left\{{\mathbf{X}}_k\right\}$ satisfying (8)–(11). Importantly, this definition does not fully specify the probability law of the process, but constrains its behavior through second-order statistical structure.

The term ${\Phi }_0{\mathbf{X}}_n$ on the right-hand side of (8) encapsulates the effect of intra-bin excitation, which induces both cross-correlations and overdispersion in the count statistics. Assuming that the spectral radius of ${\Phi }_0$ is strictly less than one, the conditional expectation and covariance matrix of ${\mathbf{X}}_n$ given $H_{n-1}$ are respectively

$$E[{\mathbf{X}}_n\mid H_{n-1}]={(I-{\Phi }_0)}^{-1}\left(\mathbf{\mu }\Delta t+\sum _{k\le n-1}{\Phi }_{n-k}{\mathbf{X}}_k\right)=:{\mathbf{\lambda }}_n^{*},$$

(12)

$$\mathrm{Var}({\mathbf{X}}_n\mid H_{n-1})={(I-{\Phi }_0)}^{-1}\mathrm{diag}\left({\mathbf{\lambda }}_n^{*}\right){(I-{\Phi }_0^{\mathsf{T}})}^{-1},$$

(13)

where I denotes the $d\times d$ identity matrix. Derivations of these expressions are provided in Appendix A.2. Due to the presence of the matrix factor ${(I-{\Phi }_0)}^{-1}$, the conditional covariance matrix is generally non-diagonal, capturing cross-correlations, and its diagonal elements exceed those of ${\mathbf{\lambda }}_n^{*}$, reflecting overdispersion. When ${\Phi }_0$ is the zero matrix, implying the absence of intra-bin excitation, the conditional variance reduces to $\mathrm{Var}({\mathbf{X}}_n\mid H_{n-1})=\mathrm{diag}\left({\mathbf{\lambda }}_n^{*}\right)$, corresponding to the conditional independence of Poisson-distributed counts.

3.3. Stationary Process

We investigate the coarse-grained Hawkes process under the assumption of stationarity and derive its second-order statistical properties. Assuming stationarity, equations (8) and (10) yield

$$\mathbf{\lambda }=E\left[{\mathbf{X}}_n\right]={\left(I-A\right)}^{-1}\mathbf{\mu }\Delta t,n\in \mathbb{Z},$$

(14)

where

$$A:=\sum _{n=0}^{\infty }{\Phi }_n$$

denotes the branching ratio matrix. Consequently, the spectral radius of A must be strictly less than one.

To derive the second-order moment of the stationary coarse-grained Hawkes process, we first establish a white noise property of the residual process.

Lemma 2.

Let $\left\{{\mathbf{X}}_n\right\}$ denote a coarse-grained Hawkes process whose branching ratio matrix has a spectral radius less than one. Then the residual process (9) forms a stationary sequence satisfying $E\left[\Delta {\mathbf{M}}_n\right]=\mathbf{0}$, $n\in \mathbb{Z}$, and

$$E\left[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_{n^{\prime }}^{\mathsf{T}}\right]={\delta }_{n{n}^{\prime }}\mathrm{diag}\left(\mathbf{\lambda }\right),n,n^{\prime }\in \mathbb{Z},$$

(15)

where ${\delta }_{n{n}^{\prime }}$ denotes the Kronecker delta.

Proof. 

See Appendix A.3.    □

As a consequence of the white noise structure of the residual process, the stationary coarse-grained Hawkes process admits a moving average representation of infinite order,

$${\mathbf{X}}_n-\mathbf{\lambda }=\sum _{k=0}^{\infty }{\Psi }_k\Delta {\mathbf{M}}_{n-k},$$

(16)

where the kernel matrices ${\Psi }_k$ are defined by

$${\Psi }_k=\sum _{j=0}^{\infty }{\left({\Phi }^{(*j)}\right)}_k,$$

with ${\Phi }^{(*j)}$ denoting the j-fold convolution of $\Phi$, and ${\left({\Phi }^{(*0)}\right)}_k={\delta }_{k0}I$. The derivation is provided in Appendix A.4. Using this representation, the autocovariance structure is obtained analogously to linear time series models.

Theorem 1.

Let $\left\{{\mathbf{X}}_n\right\}$ be a coarse-grained Hawkes process with a branching ratio matrix whose spectral radius is less than one. Then the autocovariance matrix of $\left\{{\mathbf{X}}_n\right\}$ is given by

$$R_j:=\mathrm{Cov}({\mathbf{X}}_n,{\mathbf{X}}_{n+j})=\sum _{l=0}^{\infty }{\Psi }_{l-j}\mathrm{diag}\left(\mathbf{\lambda }\right){\Psi }_l^{\mathsf{T}},j\in \mathbb{Z}.$$

Proof. 

See Appendix A.5.    □

The spectral density matrix of the stationary coarse-grained Hawkes process is obtained by taking the Fourier transform of the autocovariance sequence:

$$F_{\Delta t}^{\left(cg\right)}\left(\omega \right)={\displaystyle \frac{1}{2\pi }}\sum _{n=-\infty }^{\infty }{R}_n{e}^{-i\omega n}={\displaystyle \frac{1}{2\pi }}\widehat{\Psi }(-\omega )\mathrm{diag}\left(\mathbf{\lambda }\right){\widehat{\Psi }}^{\mathsf{T}}\left(\omega \right),$$

where $\widehat{\Psi }\left(\omega \right)={\sum }_{n=0}^{\infty }{\Psi }_n{e}^{-i\omega n}$ is the Fourier transform of the effective kernel. Alternatively, it can be expressed via the Fourier transform of the coarse-grained kernel matrix, $\widehat{\Phi }\left(\omega \right)={\sum }_{n=0}^{\infty }{\Phi }_n{e}^{-i\omega n}$, as

$$\begin{array}{c}\hfill F_{\Delta t}^{\left(cg\right)}\left(\omega \right)={\displaystyle \frac{1}{2\pi }}{(I-\widehat{\Phi }(-\omega ))}^{-1}\mathrm{diag}\left(\mathbf{\lambda }\right){(I-{\widehat{\Phi }}^{\mathsf{T}}\left(\omega \right))}^{-1},\end{array}$$

(17)

utilizing the identity $\widehat{\Psi }\left(\omega \right)={(I-\widehat{\Phi }\left(\omega \right))}^{-1}$.

In summary, under the condition that the spectral radius of the branching ratio matrix is less than one, the expected value of the process remains constant over time, and its autocovariance depends solely on the lag between observations, not their absolute positions in time. Accordingly, the coarse-grained Hawkes process is weakly stationary.

3.4. Approximation to Hawkes Process

We now establish a rigorous connection between the original Hawkes process and its coarse-grained counterpart. Specifically, we examine the second-order statistical properties of both processes under the assumption of stationarity. A summary of the statistical properties of the stationary Hawkes process is provided in Appendix B.

Since the sum of the coarse-grained kernel coincides with the integral of the excitation kernel (see Lemma 1), the branching ratio matrices, and consequently, the conditions for stationarity, are identical for both processes. Under stationarity, the expected event counts of the coarse-grained Hawkes process (14) coincide with those of the original Hawkes process (A6).

We proceed to compare the spectral density matrices of the two processes, focusing in particular on the convergence behavior of the spectral density of the coarse-grained process toward that of the original Hawkes process.

Theorem 2.

The spectral density matrix of the coarse-grained Hawkes process satisfies

$$F_{\Delta t}^{\left(cg\right)}\left(\nu \Delta t\right)=F^{\left(hw\right)}\left(\nu \right)\Delta t+O\left({\Delta t}^3\right),$$

(18)

as $\Delta t\to 0$, where $F^{\left(hw\right)}\left(\nu \right)$ denotes the spectral density matrix (A5) of the original Hawkes process.

Proof. 

See Appendix A.6.    □

Considering the binned Hawkes process, its spectral density matrix $F_{\Delta t}^{\left(hw\right)}\left(\omega \right)$ behaves as (see Appendix B)

$$F_{\Delta t}^{\left(hw\right)}\left(\nu \Delta t\right)=F^{\left(hw\right)}\left(\nu \right)\Delta t+O\left(\Delta t^3\right),$$

as $\Delta t\to 0$, which matches the expansion in (18) up to second-order terms in $\Delta t$. This leads to the following corollary.

Corollary 1.

The spectral density matrix of the coarse-grained Hawkes process approximates that of the binned Hawkes process to third-order accuracy as

$$F_{\Delta t}^{\left(cg\right)}\left(\omega \right)=F_{\Delta t}^{\left(hw\right)}\left(\omega \right)+O\left(\Delta t^3\right),$$

as $\Delta t\to 0$.

For comparative purposes, consider the binned Poisson approximation defined by equations (1)–(2). Its spectral density matrix $F_{\Delta t}^{\left(po\right)}\left(\omega \right)$ behaves as (see Appendix A.6)

$$F_{\Delta t}^{\left(po\right)}\left(\omega \right)=F_{\Delta t}^{\left(hw\right)}\left(\omega \right)+O\left(\Delta t^2\right).$$

Therefore, the coarse-grained Hawkes process yields a spectral approximation to the binned Hawkes process that is accurate to a higher order than the binned Poisson approximation.

3.5. Parameter Estimation Method

We address the problem of estimating the parameters of a Hawkes process from binned count data. Assume that we observe a sequence of binned event counts $\{{\mathbf{X}}_1,\dots ,{\mathbf{X}}_n\}$ of length n generated by a d-dimensional Hawkes process, whose excitation kernel matrix $\Phi \left(\mathbf{\theta }\right)$ is specified by a parametric form with an unknown parameter vector $\mathbf{\theta }$. We propose a parameter estimation method for $\mathbf{\theta }$ based on the coarse-grained Hawkes process framework.

To this end, we utilize the AR(∞) representation of the coarse-grained Hawkes process (see Appendix A.4),

$${\mathbf{X}}_n-\mathbf{\lambda }=\sum _{k=1}^{\infty }{\Phi }_k^{*}\left(\mathbf{\theta }\right)({\mathbf{X}}_{n-k}-\mathbf{\lambda })+\Delta {\mathbf{M}}_n^{*}\left(\mathbf{\theta }\right),$$

where ${\Phi }_k^{*}\left(\mathbf{\theta }\right)={(I-{\Phi }_0\left(\mathbf{\theta }\right))}^{-1}{\Phi }_k\left(\mathbf{\theta }\right)$, and $\Delta {\mathbf{M}}_n^{*}\left(\mathbf{\theta }\right)={(I-{\Phi }_0\left(\mathbf{\theta }\right))}^{-1}\Delta {\mathbf{M}}_n\left(\mathbf{\theta }\right)$ is a zero-mean, cross-correlated white noise sequence satisfying

$$\begin{array}{cc}\hfill E\left[\Delta {\mathbf{M}}_n^{*}\left(\mathbf{\theta }\right){\Delta {\mathbf{M}}_{n^{\prime }}^{*}\left(\mathbf{\theta }\right)}^{\mathsf{T}}\right]& ={\delta }_{n{n}^{\prime }}{(I-{\Phi }_0\left(\mathbf{\theta }\right))}^{-1}\mathrm{diag}\left(\mathbf{\lambda }\right){(I-{\Phi }_0{\left(\mathbf{\theta }\right)}^{\mathsf{T}})}^{-1}\hfill \\ \hfill & =:{\delta }_{n{n}^{\prime }}{\Lambda }^{*}\left(\mathbf{\theta }\right),n,n^{\prime }\in \mathbb{Z}.\hfill \end{array}$$

Based on this representation, we define a loss function in quadratic form, whose minimizer yields an estimator of the parameter vector:

$$\begin{array}{c}\hfill {\mathcal{L}}_n\left(\mathbf{\theta }\right)=\sum _{k=1}^n{\left({\mathbf{Y}}_k-\sum _{j=1}^{k-1}{\Phi }_j^{*}\left(\mathbf{\theta }\right){\mathbf{Y}}_{k-j}\right)}^{\mathsf{T}}{{\Lambda }^{*}}^{-1}\left(\mathbf{\theta }\right)\left({\mathbf{Y}}_k-\sum _{j=1}^{k-1}{\Phi }_j^{*}\left(\mathbf{\theta }\right){\mathbf{Y}}_{k-j}\right)+nlog\left|{\Lambda }^{*}\left(\mathbf{\theta }\right)\right|,\end{array}$$

(19)

where ${\mathbf{Y}}_k={\mathbf{X}}_k-\mathbf{\lambda }$. In practice, the unknown mean vector $\mathbf{\lambda }$ is replaced by the empirical mean $\widehat{\mathbf{\lambda }}=n^{-1}{\sum }_{k=1}^n{\mathbf{X}}_k$. The loss function then simplifies to

$$\begin{array}{c}\hfill {\mathcal{L}}_n\left(\mathbf{\theta }\right)=\sum _{k=1}^n{\left({\mathbf{Y}}_k-\sum _{j=0}^{k-1}{\Phi }_j\left(\mathbf{\theta }\right){\mathbf{Y}}_{k-j}\right)}^{\mathsf{T}}\mathrm{diag}{\left(\widehat{\mathbf{\lambda }}\right)}^{-1}\left({\mathbf{Y}}_k-\sum _{j=0}^{k-1}{\Phi }_j\left(\mathbf{\theta }\right){\mathbf{Y}}_{k-j}\right)-2nlog|I-{\Phi }_0\left(\mathbf{\theta }\right)|,\end{array}$$

(20)

where the constant term $log\mathrm{diag}\left(\widehat{\mathbf{\lambda }}\right)$ has been omitted. The first term on the right-hand side of (20) corresponds to a weighted quadratic loss, while the second term serves as a regularization component that discourages trivial solutions of the form ${\Phi }_k\left(\mathbf{\theta }\right)={\delta }_{k0}I$. The optimal parameter estimate is obtained by minimizing the loss function with respect to $\mathbf{\theta }$.

The estimation procedure is summarized as follows:

Algorithm 1: Estimation procedure for Hawkes processes

1: Given the observed bin-count sequence $\{{\mathbf{X}}_1,\dots ,{\mathbf{X}}_n\}$, compute the empirical mean $\widehat{\mathbf{\lambda }}=n^{-1}{\sum }_{k=1}^n{\mathbf{X}}_k$ and center the data as ${\mathbf{Y}}_k={\mathbf{X}}_k-\widehat{\mathbf{\lambda }}$ for $k=1,\dots ,n$. 2: Determine the parameter estimate $\widehat{\mathbf{\theta }}$ by minimizing the loss function (20). 3: Using the estimated parameter vector $\widehat{\mathbf{\theta }}$, compute the estimate of the baseline intensity as $$\widehat{\mathbf{\mu }}=(I-\widehat{A})\widehat{\mathbf{\lambda }}/\Delta ,$$ where $\widehat{A}={\sum }_{n=0}^{\infty }{\Phi }_n\left(\widehat{\mathbf{\theta }}\right)$ denotes the estimated branching ratio matrix.

4. Numerical Experiments

4.1. Assessment of Second-Order Characteristics

As established in Corollary 1, the coarse-grained Hawkes process asymptotically approximates the spectral density matrix of the binned Hawkes process as the bin size $\Delta t\to 0$. In this section, we conduct numerical investigations to evaluate the validity and robustness of this approximation for increasing values of $\Delta t$. To this end, we consider a bivariate Hawkes process ($d=2$), where each component of the excitation kernel matrix is defined as

$${\varphi }_{ij}\left(t\right)={\alpha }_{ij}{g}_{ij}\left(t\right),t\ge 0,$$

with $g_{ij}\left(t\right)$ denoting a normalized kernel satisfying ${\int }_0^{\infty }{g}_{ij}\left(t\right)dt=1$. The coefficient ${\alpha }_{ij}$ specifies the branching ratio from component j to component i, and $g_{ij}\left(t\right)$ characterizes the distribution of waiting times for event excitation. We specifically focus on a symmetric bivariate Hawkes process where the parameters satisfy ${\mu }_1={\mu }_2=\mu$, $g_{ij}\left(t\right)=g\left(t\right)$ ($1\le i,j\le 2$), and the excitation kernel matrix takes the form

$$\left(\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)=\left(\begin{array}{cc}{\alpha }^{\left(s\right)}& {\alpha }^{\left(c\right)}\\ {\alpha }^{\left(c\right)}& {\alpha }^{\left(s\right)}\end{array}\right).$$

The stationarity condition for the process holds if ${\alpha }^{\left(s\right)}+{\alpha }^{\left(c\right)}<1$. Under this constraint, the second-order statistical structure of the symmetric Hawkes process is described by the power spectral density (PSD) of each component (the diagonal elements of the spectral density matrix) and the cross-spectral density (CSD) between them (the off-diagonal elements).

We illustrate our findings using an exponential kernel defined as

$$g\left(t\right)=\beta e^{-\beta t},t\ge 0,$$

where ${\beta }^{-1}$ denotes the expected waiting time before excitation. Figure 1 presents the PSD and CSD, along with the auto- and cross-covariance functions, of the binned Hawkes process with parameters ${\alpha }^{\left(s\right)}=0.4$, ${\alpha }^{\left(c\right)}=0.3$, and $\beta =1$, for bin sizes $\Delta t=0.1$ (a), 1 (b), and 2 (c). These are depicted using blue dotted lines.

The associated coarse-grained Hawkes process is constructed using the coarse-grained kernel derived from the exponential function:

$$g_k=\left\{\begin{array}{cc}1-(1-e^{-\beta \Delta t})/\beta \Delta t,\hfill & k=0,\hfill \\ (e^{\beta \Delta t}+e^{-\beta \Delta t}-2)e^{-\beta k\Delta t}/\beta \Delta t,\hfill & k=1,2,\dots \hfill \end{array}\right.$$

In the same figures, the four second-order statistics of the coarse-grained process are represented by red lines, while those corresponding to the binned Poisson approximation are shown in green.

From these comparisons, we observe that for small bin sizes ($\Delta t=0.1<{\beta }^{-1}$), both the coarse-grained Hawkes process and the binned Poisson approximation closely reproduce the second-order behavior of the binned Hawkes process (Figure 1a). However, for $\Delta t=1$, which is comparable to the mean waiting time, the Poisson approximation significantly deteriorates (Figure 1b), whereas the coarse-grained Hawkes process continues to provide a high-fidelity approximation. Even for a larger bin size $\Delta t=2$, exceeding the characteristic time scale ${\beta }^{-1}$, the coarse-grained model remains accurate (Figure 1c).

To quantitatively assess the fidelity of the approximations, we introduce a divergence measure based on the log-determinant of the spectral density matrices. Specifically, let $F_{\Delta t}^{(••)}\left(\omega \right)$ denote the spectral density matrix of the approximating process, where `$••$’ indicates either the coarse-grained (cg) or Poisson (po) approximation. Then, the information loss relative to the binned Hawkes process $F_{\Delta t}^{\left(hw\right)}\left(\omega \right)$ is defined as

$$\Delta h^{(••)}={\displaystyle \frac{1}{4\pi }}{\int }_{-\pi }^{\pi }log{\displaystyle \frac{|F_{\Delta t}^{\left(hw\right)}\left(\omega \right)|}{|F_{\Delta t}^{(••)}\left(\omega \right)|}}d\omega ,$$

which corresponds to the gap in maximum entropy rates under spectral constraints.

Figure 2 plots the information loss for both approximations as a function of $\Delta t$, across varying values of ${\alpha }^{\left(s\right)}$ and ${\alpha }^{\left(c\right)}$. The results clearly indicate that the coarse-grained Hawkes process incurs minimal information loss across a wide range of parameter settings. In contrast, the accuracy of the binned Poisson approximation degrades with increasing $\Delta t$, with the loss exacerbated further as ${\alpha }^{\left(s\right)}$ and ${\alpha }^{\left(c\right)}$ increase.

We additionally examined the case where the excitation kernel follows a power-law distribution instead of an exponential decay. The qualitative behavior remained consistent (see Appendix C). These findings collectively demonstrate that the coarse-grained Hawkes process offers a substantially improved approximation of the second-order dynamics of the binned Hawkes process, particularly for larger bin widths.

4.2. Parameter Estimation

We now investigate the efficacy of the proposed estimation method in inferring the parameters of a Hawkes process from bin-count data. Specifically, we consider an asymmetric bivariate Hawkes process with exponential excitation kernels given by $g_{ij}\left(t\right)={\beta }_{ij}{e}^{-{\beta }_{ij}t}$ ($1\le i,j\le 2$). The parameters of the Hawkes process are set as follows:

$$\left(\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right)=\left(\begin{array}{c}1\\ 1\end{array}\right),\left(\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right)=\left(\begin{array}{cc}0.4& 0.5\\ 0.3& 0.2\end{array}\right),\left(\begin{array}{cc}{\beta }_{11}& {\beta }_{12}\\ {\beta }_{21}& {\beta }_{22}\end{array}\right)=\left(\begin{array}{cc}0.5& 0.7\\ 0.3& 1.0\end{array}\right).$$

The numerical experiments were conducted as follows. First, realizations of the Hawkes process were generated over the interval $[0,T]$ and subsequently discretized into bin-count sequences using bin size $\Delta t$. The ten model parameters ${\{{\mu }_i,{\alpha }_{ij},{\beta }_{ij}\}}_{1\le i,j\le 2}$ were then estimated from these sequences.

To assess the performance of our proposed method, we compare it against three established approaches. The first is the MC-EM algorithm proposed in [20], for which we employed the publicly available implementation [25]. The second method is maximum likelihood estimation (MLE) applied to the binned Poisson approximation. The third method involves conditional least squares estimation for the INAR(p) process, as introduced in [18]. Since the INAR(p) framework yields nonparametric estimates of the excitation kernels, we obtained the corresponding parametric kernel parameters by fitting an exponential function to the nonparametric estimates [19,20].

Figure 3 presents boxplots of the estimated values for each of the ten parameters across 500 simulated realizations of the binned Hawkes process, using $T=1000$ and $\Delta t=2$. It is apparent that both the binned Poisson MLE and the INAR(p) method produce significantly biased estimates, which is expected given that these approaches disregard excitation effects within each bin. Furthermore, we observe a substantial number of outliers in the estimates of ${\beta }_{ij}$ across all four methods, indicating that estimation of the kernel scales exhibits higher variance than estimation of the baseline intensities or branching ratios. For instance, our method yielded estimates of ${\beta }_{22}$ that deviated by a factor of 100 from the true value in 16 out of the 500 trials, highlighting the challenges in accurately estimating ${\beta }_{ij}$ when the bin size is large relative to the kernel time scale.

Figure 4 displays the root mean squared error (RMSE) of the parameter estimates across various bin sizes. For the RMSE of ${\widehat{\beta }}_{ij}$, extreme outliers with values exceeding 100 were excluded to mitigate their undue influence. Overall, the proposed method demonstrates superior performance compared to both the binned Poisson MLE and the INAR(p) approach, and achieves accuracy comparable to the MC-EM algorithm. Additionally, both the MC-EM and the proposed methods consistently maintain low RMSE values across different bin sizes, whereas the RMSEs for the binned Poisson MLE and the INAR(p) approach increase with larger bin sizes. A similar trend is observed for ${\widehat{\beta }}_{ij}$, albeit with more fluctuation due to outliers.

Figure 5 and Figure 6 illustrate, respectively, the bias and standard deviation components of the RMSE for each of the ten parameters. The proposed method yields the lowest bias among all methods, with only a few exceptions (Figure 5); meanwhile, the MC-EM algorithm attains the lowest standard deviation (Figure 6).

In conclusion, the proposed method matches the performance of the MC-EM algorithm while outperforming both the binned Poisson MLE and the INAR(p) method, particularly in terms of bias reduction. It provides robust and stable estimates for the baseline intensities ${\mu }_i$ and the branching ratios ${\alpha }_{ij}$, with estimation accuracy largely unaffected by bin size. Accurate estimation of the kernel scales ${\beta }_{ij}$ is feasible when the bin size is smaller than the characteristic kernel scale, but becomes unreliable as the bin size increases.

5. Discussion

In this study, we introduced the coarse-grained Hawkes process as an analytical approximation to the binned Hawkes process. Unlike conventional discretization techniques, the proposed framework incorporates a coarse-grained excitation kernel that systematically accounts for intra-bin excitations. Consequently, the coarse-grained Hawkes process faithfully reproduces the second-order statistical properties of the binned Hawkes process, even when the bin size exceeds the characteristic timescale of the excitation kernel. Moreover, we demonstrated that the proposed approach enables stable estimation of Hawkes process parameters from bin-count data. In particular, both the branching ratios and baseline intensities can be reliably inferred, irrespective of the temporal resolution of the bin-count sequences.

A central distinction between our approach and the Monte Carlo Expectation-Maximization (MC-EM) algorithm lies in the treatment of latent event times within the bin-count data. Whereas the MC-EM method necessitates Monte Carlo sampling from the conditional distribution of unobserved events, resulting in considerable computational burden, our method employs a parsimonious assumption that events are uniformly distributed within each bin. This assumption facilitates the analytical derivation of approximate conditional expectations. Despite its simplicity, the proposed method achieves estimation accuracy on par with that of the MC-EM algorithm while offering substantial computational advantages.

We further highlight a theoretical connection between our estimation framework and the spectral method, previously validated in the univariate setting [22]. Applying the Fourier transform, the loss function in equation (19) asymptotically approximates a spectral likelihood for large n:

$${\mathcal{L}}_n\left(\mathbf{\theta }\right)\approx \sum _{m=1}^n\left({\widehat{\mathbf{X}}}_m^{*}{F}_{\Delta t}^{\left(cg\right)}{\left({\omega }_m\right)}^{-1}{\widehat{\mathbf{X}}}_m+log\left|F_{\Delta t}^{\left(cg\right)}\left({\omega }_m\right)\right|\right),$$

where ${\omega }_m=2\pi m/n$,

$${\widehat{\mathbf{X}}}_m={\displaystyle \frac{1}{\sqrt{2\pi n}}}\sum _{k=1}^n({\mathbf{X}}_k-\mathbf{\lambda })e^{-2\pi ikm/n},$$

and ${\widehat{\mathbf{X}}}_m^{*}$ denotes the Hermitian transpose. Notably, the spectral likelihood for the multivariate binned Hawkes process can be obtained by replacing $F_{\Delta t}^{\left(cg\right)}\left({\omega }_m\right)$ with the spectral density matrix of the binned Hawkes process.

Finally, we emphasize the potential extension of the proposed framework to nonstationary time series. Owing to its formulation in the time domain, the coarse-grained Hawkes process is amenable to integration within state-space modeling paradigms [26,27], offering a promising direction for modeling nonstationary dynamics. We propose this as a compelling avenue for future research.

Appendix A.1. Proof of Lemma 1

It follows from equations (6) and (7) that

$$\begin{array}{cc}\hfill \sum _{k=0}^{\infty }{\varphi }_k& =\underset{n\to \infty }{lim}\left\{{\psi }_0+({\psi }_1-{\psi }_0)+\cdots +({\psi }_n-{\psi }_{n-1})\right\}=\underset{n\to \infty }{lim}{\psi }_n\hfill \\ \hfill & =\underset{n\to \infty }{lim}{\displaystyle \frac{1}{\Delta t}}{\int }_{n\Delta t}^{(n+1)\Delta t}\left({\int }_0^t\varphi \left(u\right)du\right)dt=\underset{n\to \infty }{lim}{\displaystyle \frac{1}{\Delta t}}{\int }_0^{\Delta t}\left({\int }_0^{n\Delta t+s}\varphi \left(u\right)du\right)ds.\hfill \end{array}$$

Since ${\int }_0^{n\Delta t+s}\varphi \left(u\right)du$ is nonnegative and monotonically nondecreasing with respect to n, we may invoke the dominated convergence theorem to interchange the limit and the integral, yielding

$$\begin{array}{c}\hfill \sum _{k=0}^{\infty }{\varphi }_k={\displaystyle \frac{1}{\Delta t}}{\int }_0^{\Delta t}\left(\underset{n\to \infty }{lim}{\int }_0^{n\Delta t+s}\varphi \left(u\right)du\right)ds={\displaystyle \frac{1}{\Delta t}}{\int }_0^{\Delta t}\left({\int }_0^{\infty }\varphi \left(u\right)du\right)ds=\alpha .\end{array}$$

Furthermore, since ${\int }_0^t\varphi \left(u\right)du$ is continuous, nonnegative, and monotonically nondecreasing in t, we obtain

$$\begin{array}{c}\hfill {\varphi }_0={\displaystyle \frac{1}{\Delta t}}{\int }_0^{\Delta t}\left({\int }_0^t\varphi \left(u\right)du\right)dt<{\displaystyle \frac{1}{\Delta t}}{\int }_0^{\Delta t}\left({\int }_0^{\infty }\varphi \left(u\right)du\right)dt=\alpha .\end{array}$$

Appendix A.2. Derivation of (uid28) and (uid29)

From equations (8) and (10), we have

$$\begin{array}{cc}\hfill E[{\mathbf{X}}_n\mid H_{n-1}]& =E\left[\mathbf{\mu }\Delta t+\sum _{k\le n}{\Phi }_{n-k}{\mathbf{X}}_k\mid H_{n-1}\right]\hfill \\ \hfill & =\mathbf{\mu }\Delta t+\sum _{k\le n-1}{\Phi }_{n-k}{\mathbf{X}}_k+{\Phi }_{0}E[{\mathbf{X}}_n\mid H_{n-1}].\hfill \end{array}$$

$$\therefore (I-{\Phi }_0)E[{\mathbf{X}}_n\mid H_{n-1}]=\mathbf{\mu }\Delta t+\sum _{k\le n-1}{\Phi }_{n-k}{\mathbf{X}}_k.$$

Since the spectral radius of ${\Phi }_0$ is strictly less than one, the matrix $(I-{\Phi }_0)$ is invertible, and all entries of $E[{\mathbf{X}}_n\mid H_{n-1}]$ are nonnegative. Hence, the conditional mean vector of ${\mathbf{X}}_n$ is given by expression (12).

To derive (13), observe that

$$\begin{array}{c}\hfill {\mathbf{X}}_n-{\mathbf{\lambda }}_n^{*}={(I-{\Phi }_0)}^{-1}\left((I-{\Phi }_0){\mathbf{X}}_n-\mathbf{\mu }\Delta t-\sum _{k\le n-1}{\Phi }_{n-k}{\mathbf{X}}_k\right)={(I-{\Phi }_0)}^{-1}\Delta {\mathbf{M}}_n.\end{array}$$

Therefore,

$$\begin{array}{cc}\hfill \mathrm{Var}({\mathbf{X}}_n\mid H_{n-1})& =E\left[({\mathbf{X}}_n-{\mathbf{\lambda }}_n^{*}){({\mathbf{X}}_n-{\mathbf{\lambda }}_n^{*})}^{\mathsf{T}}\mid H_{n-1}\right]\hfill \\ \hfill & ={(I-{\Phi }_0)}^{-1}E\left[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_n^{\mathsf{T}}\mid H_{n-1}\right]{(I-{\Phi }_0^{\mathsf{T}})}^{-1}\hfill \\ \hfill & ={(I-{\Phi }_0)}^{-1}\mathrm{diag}\left(E[{\mathbf{\lambda }}_n\mid H_{n-1}]\right){(I-{\Phi }_0^{\mathsf{T}})}^{-1},\hfill \end{array}$$

where the last equality follows from equation (11). Finally, since $E[{\mathbf{\lambda }}_n\mid H_{n-1}]={\mathbf{\lambda }}_n^{*}$ by (10), we obtain equation (13).

Appendix A.3. Proof of Lemma 2

From equation (10), it follows that

$$E\left[\Delta {\mathbf{M}}_n\right]=E\left[E[\Delta {\mathbf{M}}_n\mid H_{n-1}]\right]=\mathbf{0}.$$

To derive equation (15), observe that for $n<n^{\prime }$ (and, by symmetry, for $n\ne n^{\prime }$),

$$E\left[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_{n^{\prime }}^{\mathsf{T}}\right]=E\left[\Delta {\mathbf{M}}_{n}E[\Delta {\mathbf{M}}_{n^{\prime }}^{\mathsf{T}}\mid H_{n^{\prime }-1}]\right]=\mathbf{0}.$$

Moreover, applying equation (11), we obtain

$$\begin{array}{c}\hfill E\left[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_n^{\mathsf{T}}\right]=E\left[E[\Delta {\mathbf{M}}_n\Delta {\mathbf{M}}_n^{\mathsf{T}}\mid H_{n-1}]\right]=E\left[\mathrm{diag}\left(E[{\mathbf{\lambda }}_n\mid H_{n-1}]\right)\right]=\mathrm{diag}\left(\mathbf{\lambda }\right),\end{array}$$

which completes the proof.

Appendix A.4. AR(∞) and MA(∞) representations

Using equations (9) and (14), the rate equation (8) can be reformulated as

$${\mathbf{X}}_n-\mathbf{\lambda }=\sum _{k=0}^{\infty }{\Phi }_k({\mathbf{X}}_{n-k}-\mathbf{\lambda })+\Delta {\mathbf{M}}_n,$$

(A1)

$$\therefore (I-{\Phi }_0)({\mathbf{X}}_n-\mathbf{\lambda })=\sum _{k=1}^{\infty }{\Phi }_k({\mathbf{X}}_{n-k}-\mathbf{\lambda })+\Delta {\mathbf{M}}_n.$$

Under the stationarity condition, the inverse ${(I-{\Phi }_0)}^{-1}$ exists, yielding the autoregressive representation

$${\mathbf{X}}_n-\mathbf{\lambda }=\sum _{k=1}^{\infty }{\Phi }_k^{*}({\mathbf{X}}_{n-k}-\mathbf{\lambda })+\Delta {\mathbf{M}}_n^{*},$$

where ${\Phi }_k^{*}={(I-{\Phi }_0)}^{-1}{\Phi }_k$ and $\Delta {\mathbf{M}}_n^{*}={(I-{\Phi }_0)}^{-1}\Delta {\mathbf{M}}_n$.

To obtain a moving average representation, we apply the formal z-transform to both sides of equation (A1):

$$\widehat{\mathbf{Y}}\left(z\right)=\widehat{\Phi }\left(z\right)\widehat{\mathbf{Y}}\left(z\right)+\Delta \widehat{\mathbf{M}}\left(z\right),$$

where $\widehat{\mathbf{Y}}\left(z\right)={\sum }_{n=-\infty }^{\infty }({\mathbf{X}}_n-\mathbf{\lambda })z^{-n}$, $\widehat{\Phi }\left(z\right)={\sum }_{n=-\infty }^{\infty }{\Phi }_n{z}^{-n}$ and $\Delta \widehat{\mathbf{M}}\left(z\right)={\sum }_{n=-\infty }^{\infty }\Delta {\mathbf{M}}_n{z}^{-n}$. Solving for $\widehat{\mathbf{Y}}\left(z\right)$ yields

$$\widehat{\mathbf{Y}}\left(z\right)=\widehat{\Psi }\left(z\right)\Delta \widehat{\mathbf{M}}\left(z\right),$$

where the transfer function is defined by

$$\widehat{\Psi }\left(z\right)={(I-\widehat{\Phi }\left(z\right))}^{-1}=\sum _{j=0}^{\infty }{\widehat{\Phi }}^j\left(z\right).$$

Applying the inverse z-transform then gives the MA(∞) process:

$${\mathbf{X}}_n-\mathbf{\lambda }=\sum _{k=0}^{\infty }{\Psi }_k\Delta {\mathbf{M}}_{n-k},$$

where ${\Psi }_k={\sum }_{j=0}^{\infty }{\left({\Phi }^{(*j)}\right)}_k$.

Appendix A.5. Proof of Theorem 1

Using equations (15) and (16), the autocovariance matrix of $\left\{{\mathbf{X}}_n\right\}$ is given by

$$\begin{array}{cc}\hfill \mathrm{Cov}({\mathbf{X}}_n,{\mathbf{X}}_{n+j})& =E\left[({\mathbf{X}}_n-\mathbf{\lambda }){({\mathbf{X}}_{n+j}-\mathbf{\lambda })}^{\mathsf{T}}\right]\hfill \\ \hfill & =E\left[\left(\sum _{k=0}^{\infty }{\Psi }_k\Delta {\mathbf{M}}_{n-k}\right){\left(\sum _{l=0}^{\infty }{\Psi }_l\Delta {\mathbf{M}}_{n+j-l}\right)}^{\mathsf{T}}\right]\hfill \\ \hfill & =\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }{\Psi }_{k}E\left[\Delta {\mathbf{M}}_{n-k}\Delta {\mathbf{M}}_{n+j-l}^{\mathsf{T}}\right]{\Psi }_l^{\mathsf{T}}\hfill \\ \hfill & =\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }{\Psi }_k\mathrm{diag}\left(\mathbf{\lambda }\right){\delta }_{k,l-j}{\Psi }_l^{\mathsf{T}}\hfill \\ \hfill & =\sum _{l=0}^{\infty }{\Psi }_{l-j}\mathrm{diag}\left(\mathbf{\lambda }\right){\Psi }_l^{\mathsf{T}},\hfill \end{array}$$

which completes the proof.

Appendix A.6. Proof of Theorem 2

We provide a proof of Theorem 2 in the univariate setting. The extension to the multivariate case follows by applying the same reasoning component-wise to each element of the spectral density matrix. To begin, we establish the following approximation result for the Fourier transforms of the excitation kernel and its coarse-grained counterpart.

Lemma A1.

Let $\tilde{\varphi }\left(\nu \right)={\int }_0^{\infty }\varphi \left(t\right)e^{-i\nu t}dt$ denote the continuous-time Fourier transform of the excitation kernel, and let $\widehat{\varphi }\left(\omega \right)={\sum }_{k=0}^{\infty }{\varphi }_k{e}^{-i\omega k}$ denote the discrete-time Fourier transform of the corresponding coarse-grained kernel. Then, in the limit as $\Delta t\to 0$, the following approximation holds:

$$\widehat{\varphi }\left(\nu \Delta t\right)=\tilde{\varphi }\left(\nu \right)+O\left(\Delta t^2\right).$$

(A2)

Proof. 

From (6) and (7), we obtain

$$\begin{array}{cc}\hfill \widehat{\varphi }\left(\nu \Delta t\right)& =\sum _{k=0}^{\infty }{\varphi }_k{e}^{-i\nu k\Delta t}\hfill \\ \hfill & ={\psi }_0+({\psi }_1-{\psi }_0)e^{-i\nu \Delta t}+({\psi }_2-{\psi }_1)e^{-2i\nu \Delta t}+\cdots \hfill \\ \hfill & =\underset{n\to \infty }{lim}\left[\sum _{k=0}^{n-1}{\psi }_k{e}^{-i\nu k\Delta t}(1-e^{-i\nu \Delta t})+{\psi }_n{e}^{-i\nu n\Delta t}\right]\hfill \\ \hfill & =\underset{n\to \infty }{lim}\left[\sum _{k=0}^{n-1}{\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu k\Delta t}dt·{\displaystyle \frac{1-e^{-i\nu \Delta t}}{\Delta t}}+{\psi }_n{e}^{-i\nu n\Delta t}\right],\hfill \end{array}$$

(A3)

where $\phi \left(t\right)={\int }_0^t\varphi \left(u\right)du$. To evaluate the integral, consider

$$\begin{array}{cc}\hfill {\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu t}dt& ={\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)\left[e^{-i\nu k\Delta t}-i\nu e^{-i\nu k\Delta t}(t-k\Delta t)+O\left(\Delta t^2\right)\right]dt\hfill \\ \hfill & ={\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu k\Delta t}dt-Q_k+O\left(\Delta t^3\right),\hfill \end{array}$$

(A4)

where

$$\begin{array}{c}\hfill Q_k=i\nu e^{-i\nu k\Delta t}{\int }_0^{\Delta t}\phi (t+k\Delta t)tdt.\end{array}$$

By Taylor’s theorem, for each $t\in [0,\Delta t]$, there exists $c\in (k\Delta t,(k+1\left)\Delta t\right)$ such that

$$\phi (t+k\Delta t)=\phi \left(k\Delta t\right)+\dot{\phi }\left(c\right)t=\phi \left(k\Delta t\right)+\varphi \left(c\right)t.$$

Note that c generally depends on t. Then $Q_k$ is approximated as

$$\begin{array}{c}\hfill Q_k=i\nu e^{-i\nu k\Delta t}{\int }_0^{\Delta t}\left[\phi \left(k\Delta t\right)+\varphi \left(c\right)t\right]tdt=i\nu e^{-i\nu k\Delta t}·{\displaystyle \frac{\phi \left(k\Delta t\right)}{2}}\Delta t^2+O\left(\Delta t^3\right).\end{array}$$

Substituting into (A4), we obtain

$$\begin{array}{c}\hfill {\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu k\Delta t}dt={\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu t}dt+i\nu e^{-i\nu k\Delta t}·{\displaystyle \frac{\phi \left(k\Delta t\right)}{2}}\Delta t^2+O\left({\Delta }^3\right).\end{array}$$

Therefore, for $n\gg 1$ and $\Delta t\sim n^{-1}$, it follows that

$$\begin{array}{cc}\hfill \sum _{k=0}^{n-1}{\int }_{k\Delta t}^{(k+1)\Delta t}\phi \left(t\right)e^{-i\nu k\Delta t}dt& ={\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+{\displaystyle \frac{i\nu \Delta t}{2}}\sum _{k=0}^{n-1}\phi \left(k\Delta t\right)e^{-i\nu k\Delta t}\Delta t+O\left(\Delta t^2\right)\hfill \\ \hfill & ={\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+{\displaystyle \frac{i\nu \Delta t}{2}}{\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+O\left(\Delta t^2\right).\hfill \end{array}$$

Substituting back into (A3) yields

$$\begin{array}{cc}\hfill \widehat{\varphi }\left(\nu \Delta t\right)& =\underset{n\to \infty }{lim}[\left({\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+{\displaystyle \frac{i\nu \Delta t}{2}}{\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+O\left(\Delta t^2\right)\right)\hfill \\ \hfill & \times \left(i\nu +{\displaystyle \frac{{\nu }^2\Delta t}{2}}+O\left(\Delta t^2\right)\right)+{\psi }_n{e}^{-i\nu n\Delta t}]\hfill \\ \hfill & =\underset{n\to \infty }{lim}\left(i\nu {\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt+{\psi }_n{e}^{-i\nu n\Delta t}\right)+O\left(\Delta t^2\right).\hfill \end{array}$$

The first term on the right-hand side can be further simplified as

$$\begin{array}{cc}\hfill i\nu {\int }_0^{n\Delta t}\phi \left(t\right)e^{-i\nu t}dt& =-{\left[\phi \left(t\right)e^{-i\nu t}\right]}_0^{n\Delta t}+{\int }_0^{n\Delta t}\dot{\phi }\left(t\right)e^{-i\nu t}dt\hfill \\ \hfill & =-\phi \left(n\Delta t\right)e^{-i\nu n\Delta t}+{\int }_0^{n\Delta t}\varphi \left(t\right)e^{-i\nu t}dt.\hfill \end{array}$$

Consequently,

$$\begin{array}{cc}\hfill \widehat{\varphi }\left(\nu \Delta t\right)& =\underset{n\to \infty }{lim}\left({\int }_0^{n\Delta t}\varphi \left(t\right)e^{-i\nu t}dt-\phi \left(n\Delta t\right)e^{-i\nu n\Delta t}+{\psi }_n{e}^{-i\nu n\Delta t}\right)+O\left(\Delta t^2\right)\hfill \\ \hfill & ={\int }_0^{\infty }\varphi \left(t\right)e^{-i\nu t}dt+O\left(\Delta t^2\right),\hfill \end{array}$$

where we have used the identity ${lim}_{n\to \infty }{\psi }_n=\phi (\infty )$. □

Now we complete a proof of Theorem 2. According to equation (17), the spectral density of the univariate coarse-grained Hawkes process is expressed as

$$f_{\Delta t}^{\left(cg\right)}\left(\omega \right)={\displaystyle \frac{\lambda \Delta t}{2\pi }}{\displaystyle \frac{1}{|1-\widehat{\varphi }{\left(\omega \right)|}^2}}.$$

Substituting equation (A2) into the expression above yields

$$\begin{array}{cc}\hfill f_{\Delta t}^{\left(cg\right)}\left(\nu \Delta t\right)& ={\displaystyle \frac{\lambda \Delta t}{2\pi }}{\displaystyle \frac{1}{|1-\tilde{\varphi }\left(\nu \right)+O\left(\Delta t^2\right){|}^2}}\hfill \\ \hfill & ={\displaystyle \frac{\lambda \Delta t}{2\pi }}\left({\displaystyle \frac{1}{|1-\tilde{\varphi }{\left(\nu \right)|}^2}}+O\left(\Delta t^2\right)\right)\hfill \\ \hfill & =f^{\left(hw\right)}\left(\nu \right)\Delta t+O\left(\Delta t^3\right),\hfill \end{array}$$

thereby completing the proof of Theorem 2 in the univariate setting.

For comparison, consider the binned Poisson approximation introduced in equations (1)–(2), wherein the discretized excitation kernel is defined by ${\varphi }_k=\varphi \left(k\Delta t\right)\Delta t$ for $k=1,2,\dots$. The discrete-time Fourier transform of this kernel is given by

$$\begin{array}{c}\hfill \widehat{\varphi }\left(\nu \Delta t\right)=\sum _{k=1}^{\infty }\varphi \left(k\Delta t\right)e^{-i\nu k\Delta t}\Delta t=\tilde{\varphi }\left(\nu \right)+O\left(\Delta t\right),\end{array}$$

as $\Delta \to 0$. Observe that the approximation order is lower than that of equation (A2). Consequently, the spectral density of the binned Poisson approximation is asymptotically given by

$$f_{\Delta t}^{\left(po\right)}\left(\nu \Delta t\right)=f^{\left(hw\right)}\left(\nu \right)\Delta t+O\left(\Delta t^2\right).$$