# Multivariate Representations of Univariate Marked Hawkes Processes

Louis Davis1,2 , Conor Kresin1, Boris Baeumer1, Ting Wang1
1Department of Mathematics and Statistics, University of Otago,
P.O. Box 56, Dunedin 9054, New Zealand
2E-mail<EMAIL_ADDRESS>
###### Abstract

Univariate marked Hawkes processes are used to model a range of real-world
phenomena including earthquake aftershock sequences, contagious disease
spread, content diffusion on social media platforms, and order book dynamics.
This paper illustrates a fundamental connection between univariate marked
Hawkes processes and multivariate Hawkes processes. Exploiting this connection
renders a framework that can be built upon for expressive and flexible
inference on diverse data. Specifically, multivariate unmarked Hawkes
representations are introduced as a tool to parameterize univariate marked
Hawkes processes. We show that such multivariate representations can
asymptotically approximate a large class of univariate marked Hawkes
processes, are stationary given the approximated process is stationary, and
that resultant conditional intensity parameters are identifiable. A simulation
study demonstrates the efficacy of this approach, and provides heuristic
bounds for error induced by the relatively larger parameter space of
multivariate Hawkes processes.

Keywords: Marked point process; Multivariate Hawkes Process; Non-separable
conditional intensity

##  1 Introduction

The Hawkes point process model [18] is widely used to model branching
processes of interest ranging from contagious disease spread [39, 45] and
retaliatory crime [28] to earthquake aftershock sequences [31], neural spike
train data [6], popularity of content on social media [40], and high-frequency
trading order book characteristics [3]. Hawkes models and their variants are
often employed to capture temporal clustering between points, roughly serving
as an analog to autoregressive models for time series data that does not
suffer from the modifiable areal unit problem [12, 17]. Hawkes models can
accommodate marked point process data, which includes spatiotemporal data and
event stream data as two categories of particular interest [12, 15], as well
as multivariate data wherein different types of points can mutually excite or
inhibit [9, 24].

First order conditional intensity specifications of univariate marked Hawkes
processes often assume mark-separability, an analytically and computationally
necessary model mis-specification. For instance, the ground intensity and
magnitude distribution of earthquake aftershocks are modelled as separable
using the Epidemic Type Aftershock Sequence (ETAS) model [31] despite the fact
that the ground intensity and magnitude distribution are postulated to be non-
separable [14, 42, 50]. Although non-separable marked Hawkes processes provide
a very general modelling framework, they are rarely implemented, as model
specification is frequently challenging in practice [53].

Parameters for conditional intensity functions governing point processes are
often estimated via the method of maximum likelihood (MLE), as the resultant
estimates are consistent and asymptotically normal [29]. The properties of MLE
estimates of univariate Hawkes processes are discussed in detail in [36], and
are often conjectured to hold for marked univariate and multivariate Hawkes
processes [7]. Reliance on MLE estimates increases the practical difficulty of
fitting marked univariate Hawkes models as convergence to parameter estimates
is slow [47] and the computational expense of likelihood calculation is high
[37]. Less computationally costly estimation methods such as expectation
maximization (EM) [47], spectral estimation [10], and Takacs Fiksel estimation
[22] are available, but are relatively unstudied and infrequently implemented
for estimation of marked Hawkes processes.

Precedent work demonstrating that a non-separable univariate marked Hawkes
process can be expressed as a multivariate marked Hawkes process has shown
promising results [14]. However, the difficulty of computation and model
specification of fitting multivariate marked Hawkes processes is exacerbated
relative to univariate marked Hawkes processes. In the context of simulation
of multivariate Hawkes processes via thinning, the foundational connection
between multivariate and univariate Hawkes processes is well developed, see
for instance Proposition 1 of [30]. Further, the more general theoretical
framework of predictable projection of multivariate point processes is
detailed in [20].

We propose partitioning the support of the mark distribution of a univariate
marked Hawkes process and then modelling the process as an unmarked
multivariate Hawkes process such that each component intensity corresponds to
a distinct mark range. Our approach allows for both seperable and non-
separable marked Hawkes processes to be captured flexibly via cross excitation
in the multivariate Hawkes representation. The granularity of cross-excitation
with respect to the mark distribution can be specified via the selected
partitioning of the mark space, allowing for the true mark distribution to be
unknown and free of modelling assumptions. Further, our approach retains the
benefits of modelling the univariate marked process parametrically as the
estimated parameters are interpretable.

Two existing methods utilise multivariate Hawkes representations of non-
separable univariate marked Hawkes processes. The first existing approach is a
class of Hawkes-like models, here termed Neural Hawkes. Neural Hawkes models
are neural networks built on architectures ranging from continuous long short-
term memory (cLSTM) [27] to transformer [55]. Such models are commonly applied
in the case of event stream data; such data is naturally modelled as a non-
separable univariate marked Hawkes process with a discrete mark space (wherein
each mark is a type of point). Neural Hawkes models are flexible unmarked
multivariate models that can parameterize cross excitement, model
nonstationarity, nonlinearity, and even negative background rates. However,
such models do not yield interpretable parameters and the finite sample and
asymptotic statistical properties of Neural Hawkes estimates are yet to be
developed [41, 46]. In contrast, our method relies on the well-trodden theory
surrounding multivariate Hawkes processes [24].

A second approach, closer to our proposed method, partitions the mark space of
a (possibly non-separable) univariate marked Hawkes process whereupon the
second order characteristics of the resultant multivariate Hawkes
representation are shown to be the solution to a system of Fredholm integrals
estimable via the Wiener-Hopf method [4]. This methodology has been applied to
order book data [26], and allows for both exciting and inhibitory effects non-
parametrically. Unlike our approach, this methodology relies on completely
nonparametric estimation of the kernel matrix, and like the cLSTM, results in
non-interpretable estimates. Further, the multivariate Hawkes representation
used in our method is shown to approximate a univariate marked Hawkes process
arbitrarily well; both above-mentioned existing methods fail to provide any
such guarantee.

###  1.1 Our contribution

We show that for any given realization, there exists a multivariate unmarked
Hawkes process representation that can approximate the conditional intensity
function of a univariate marked Hawkes process so that the $L^{1}$ difference
between the two decreases as the number of components increases. Furthermore,
we demonstrate that such a multivariate representation is stationary given
stationarity of the targeted univariate marked Hawkes process. Finally, we
demonstrate that the parameters governing the multivariate Hawkes
representation are identifiable. This method excels relative to other existing
options with respect to the cases where (1) the true mark distribution is
unknown or non-stationary, and (2) if the mark distribution and ground
intensity of the univariate marked process are non-separable or require
parametric model specifications that are computationally intractable. Cross
excitement parameterized by a multivariate Hawkes process allows for the
ground intensity and mark distribution interaction present in a marked
univariate Hawkes process to be flexibly modelled in a computationally
tractable setting. These results are demonstrated in a very general setting
wherein we only assume the existence of a positive mark density, continuity of
mark-variable productivity function, and integrability of the positive kernel
function.

The remainder of the paper is structured as follows: Section 2 contains
necessary notation and preliminaries, Section 3 contains the above-described
results, Section 4 contains a simulation study validating the results
described in Section 3, and Section 5 concludes and briefly discusses future
work.

##  2 Introduction to Point Processes and Hawkes Processes

A point process $N$ is a measurable mapping from a filtered probability space
$(\Omega,\mathcal{H},\mathcal{P})$ onto $\mathcal{N}$, the set of
$\mathbb{Z}^{+}$-valued random measures (counting measures) on a complete
separable metric space (CSMS) [12]. We will restrict our attention to point
processes that are boundedly finite, i.e. processes having only a finite
number of points inside any bounded set.

Such a point process is assumed to be adapted to the filtration
$\\{\mathcal{H}_{t}\\}_{t\geq 0}$ containing all information on the process
$N$ at all locations and all times up to and including time $t$. A process is
$\mathcal{H}$-predictable if it is adapted to the filtration generated by the
left continuous processes $\mathcal{H}_{t}$. Intuitively, $\mathcal{H}_{t}$
represents the history of a process up to, but not including time $t$. A
rigorous definition of $\mathcal{H}_{t}$ can be found in [13]. For a Borel set
$X$, we say that $N(X)$ is the number of the points in the set $X$, and
$N(T)\vcentcolon=N\left([0,T)\right)$ is the number of points that have
arrived by time $T$. Assuming it exists, the $\mathcal{H}$-conditional
intensity $\lambda$ of $N$ is an integrable, non-negative,
$\mathcal{H}$-predictable process such that

$\lambda(t,M|\mathcal{H}_{t})=\lim_{h,\delta\downarrow
0}\frac{\mathbb{E}[N\left([t,t+h)\times\mathbb{B}(M,\delta)\right)|\mathcal{H}_{t}]}{h\mu\left(\mathbb{B}(M,\delta)\right)}.$
(1)

where $\mathbb{B}(M,\delta)$ is a ball centered at $M$ with radius $\delta$,
$\mu$ is the measure on $\mathcal{M}$ and $\mathcal{H}_{t}$ represents the
history of the process $N$ up to but not including time $t$.

A point process is simple if with probability one, all the points are
distinct. Since the conditional intensity $\lambda$ uniquely determines the
finite-dimensional distributions of any simple point process (Proposition
7.2.IV of [12]), one typically models a simple point process by specifying a
model for $\lambda$. A point process is stationary if for all bounded Borel
subsets of the real line $\mathcal{B}_{r}$ and times $t\in\mathbb{R}$,
$\\{N(\mathcal{B}_{r}+t)\\}_{r=1}^{\infty}\perp t$, which is to say the joint
distribution of the specified model has a structure which is invariant over
shifts in time.

A multivariate Hawkes process [18] is a simple point process
$\\{N_{i}(t)\\}_{i=1}^{k}$ or equivalently $\bm{N}(t)$, which can be
represented as a branching process. The process is defined as the sum of the
$K$ component intensities, each expressed as a Stieltjes integral

$\lambda^{(i)}(t|\mathcal{H}_{t})=\lambda_{0i}+\sum_{j=1}^{K}\int_{0}^{t}g_{ij}(t-u)dN_{j}(u).$
(2)

When $K=1$, such a process is referred to as a univariate Hawkes process. The
properties of general multivariate processes [8] and nonlinear Hawkes
processes [9] are well developed, but the scope of this paper is restricted to
linear Hawkes processes. Such processes can be viewed as Poisson cluster
processes made up of immigrant events of type $i$, arriving as a Poisson
process with rate $\lambda_{0i}$ and offspring events corresponding to each
immigrant event given a history (e.g. [1]). Offspring events are those excited
by previous arrivals in the history of the process, and we say that event
$t_{i}$ is a first generation offspring of event $t_{j}$, if $t_{j}$ is the
parent of $t_{i}$. Hence, Hawkes processes are stationary if
$\rho\left(\int_{0}^{\infty}\bm{g}(u)d\bm{N}(u)\right)<1$ where $\rho(\cdot)$
denotes the spectral radius. In this case, the mean cluster size is finite
(Lemma 6.3.II of [12]).

The conditional intensity function of a marked univariate Hawkes process takes
the parametric form

$\lambda(t,M|\mathcal{H}_{t})=\lambda_{0}(M)+\int_{0}^{t}g(t-u,M)dN(u).$ (3)

If the factorization $g(t-u,M)=f(M)g^{\prime}(t-u)$, $\lambda_{0}(M)=\Lambda
f(M)$ exists, where $f(\cdot)$ is the mark density, $g^{\prime}(\cdot)$ is the
kernel function, and $\Lambda$ is the Poisson immigrant rate, then the process
characterized in Equation (3) is mark separable.

Within this article we are interested in the properties of a multivariate
Hawkes process on the space $[0,T]\times\mathcal{M}$ where $\mathcal{M}$ is a
compact CSMS equipped with metric $d$ and measure $\mu$ such that
$\mu(\mathcal{M})<\infty$. When comparing our multivariate process to a marked
univariate process we use the same history for each process. Therefore, our
convergence theorems can be viewed as pointwise deterministic. It is likely
that our process converges to marked univariate Hawkes processes in the
Skorokhod topology, however this is outside of the scope of the paper and not
pertinent for practical purposes.

We now introduce our multivariate Hawkes representation which will be used to
approximate a non-separable univariate marked Hawkes process. Specifically, we
define the conditional intensity function, dependent on the parameter vector
$\theta$ to be

$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\sum_{i=1}^{K}\chi_{A_{i}}(M)\left(\lambda_{0i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j};\beta_{ij})\right)$
(4)

where $\chi_{A}(M)$ is the indicator function for the measurable set
$A\subset\mathcal{M}$, $\lambda_{0i}$ is the background rate parameter,
$t_{\ell,j}$ is the $\ell^{th}$ event of type $j$ (i.e. $M_{\ell,j}\in A_{j}$)
and $g_{ij}(t;\beta_{ij})$ is a positive density dependent on only one
parameter $\beta_{ij}$, which we now shorten to $g_{ij}(t)$ to ease notation.
For any fixed $K$ the parameter vector is
$\theta=\\{\\{\lambda_{0i}\\},\\{\alpha_{ij}\\},\\{\beta_{ij}\\}\\}_{{1\leq
i,j\leq K}}$, hence there are $2K^{2}+K$ parameters, thus
$\theta\in\Theta\subset\mathbb{R}^{K(2K+1)}_{+}$.

Such a representation is constructed on a set of measurable sets
$\left\\{A_{j}\right\\}_{1\leq j\leq K}$ that partition $\mathcal{M}$ such
that

  * P1.

$K\in\mathbb{N}$ is finite,

  * P2.

$A_{i}\cap A_{j}=\emptyset$ if $i\neq j$,

  * P3.

$\bigcup_{i=1}^{K}A_{i}=\mathcal{M}$, and

  * P4.

$\mu(A_{i})>0$ for every $i=1,2,\dots K$.

The purpose of properties P1-P3 is to construct a finite partition of
$\mathcal{M}$. P4 ensures that the probability an event has a mark taking a
value in $A_{i}$ is positive for every set $A_{i}$ and so that the density has
a finite mean value on the set $A_{i}$. For example, $\mathcal{M}$ could be a
compact subset of $\mathbb{R}^{d}$, in which case we would use the usual
Euclidean metric and Lebesgue measure, and the point process would then be
spatial-temporal. Spatial-temporal Hawkes processes have many applications
such as earthquakes [32, 34], crime [28] or terrorism [21]. Alternatively, if
$\mathcal{M}\subset\mathbb{Z}^{d}$ one should use the discrete metric and
counting measure. Hawkes processes with integer marks have been previously
applied to social media event stream data [38], graphical representations of
multivariate Hawkes processes [15] and a wide variety of financial data
applications [16].

##  3 Results

In this section we present the main results of our article in full detail.
Specifically, we demonstrate the existence of a multivariate representation
that approximates a univariate marked Hawkes process in the $L^{1}$ sense.
Following this, under slightly stricter assumptions, we also demonstrate that
the multivariate representation is stationary if the target process is also
stationary. We finally demonstrate that the parameters of the multivariate
representation are identifiable.

###  3.1 The special case of Mark-Separable Point Processes

In this subsection we first demonstrate that there exists a multivariate
representation of a compound Poisson process. We demonstrate that the MLE of
the multivariate representation, without excitation, are exactly those of a
histogram estimator. This specific case is included to build intuition as to
what our multivariate representation does. To do so, we require Proposition
14.3.II (c) of [13].

######  Proposition 3.1 (Proposition 14.3.II (c) of [13]).

Given a predictable conditional intensity $\lambda(t,M|\mathcal{H}_{t})$, the
associated ground process has predictable conditional intensity

$\lambda_{g}(t|\mathcal{H}_{t})=\int_{\mathcal{M}}\lambda(t,M|\mathcal{H}_{t})d\mu(M)$
(5)

and the associated mark distribution with density

$f(M|\mathcal{H}_{t})=\frac{\lambda(t,M|\mathcal{H}_{t})}{\lambda_{g}(t|\mathcal{H}_{t})}.$
(6)

We use Proposition 3.1 in the following example to demonstrate the use of our
multivariate representation.

######  Example 1.

Consider a compound Poisson process with conditional intensity
$\Lambda_{0}f(M)$ where $\Lambda_{0}>0$ is the arrival rate of events and
$f(M)$ is the density for the independent and identically distributed marks
taking values in a compact subset of $\mathbb{R}^{d}$. Suppose
$\left\\{A_{j}\right\\}_{1\leq j\leq K}$ satisfy properties P1-P4 and that
$\alpha_{ij}=0$ for every $i,j$. In this case, the MLE of the multivariate
representation (Equation (4)) are the same as the estimates given by a
histogram estimator.

To demonstrate this, let $\\{(t_{i},M_{i})\\}_{1\leq i\leq
N(T)}\subset[0,T]\times\mathcal{M}$ be a realisation of a compound Poisson
process. By Proposition 7.3.III of [12] the log-likelihood for the
multivariate representation of the compound Poisson process is

$\displaystyle\mathcal{L}(\bm{\theta})=$
$\displaystyle\sum_{j=1}^{N(T)}\log\left(\sum_{i=1}^{K}\chi_{A_{i}}(M_{j})\lambda_{0i}\right)-\int_{\mathcal{M}}\int_{0}^{T}\sum_{i=1}^{K}\chi_{A_{i}}(M)\lambda_{0i}dtd\mu(M)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{N(T)}\log(\lambda_{0k_{j}})-\sum_{i=1}^{K}\mu(A_{i})\lambda_{0i}T=\sum_{i=1}^{K}\left(N_{i}(T)\log(\lambda_{0i})-\lambda_{0i}\mu([0,T]\times
A_{i})\right)$ (7)

where $k_{j}$ is the type of event $j$ (i.e. $M_{j}\in A_{k_{j}}$), $N_{i}(T)$
is the number of events that have arrived by time $T$ such that $M\in A_{i}$
and $\mu$ is the Lebesgue measure on
$[0,T]\times\mathcal{M}\subset\mathbb{R}^{+}\times\mathbb{R}^{d}$. Since each
summand on the right hand side of Equation (7) is independent, $\lambda_{0i}$
can be estimated by maximising each summand separately. Define

$f(\lambda_{0i})=N_{i}(T)\log(\lambda_{0i})-\mu([0,T]\times
A_{i})\lambda_{0i}.$

It directly follows that

$f^{\prime}(\lambda_{0i})=\frac{N_{i}(T)}{\lambda_{0i}}-\mu([0,T]\times
A_{i})=0\implies\hat{\lambda}_{0i}=\frac{N_{i}(T)}{\mu([0,T]\times
A_{i})}=\frac{N_{i}(T)}{T\mu(A_{i})}$

which is clearly a maximum by taking the second derivative of $f$.
Furthermore, $\hat{\lambda}_{0i}$ is exactly the MLE for a Poisson process
observed on $[0,T]\times A_{i}$, therefore inheriting consistency and
asymptotic normality, e.g. Example 1 of [29].

We now examine the induced ground process for the estimated multivariate
representation using Proposition 3.1. By Equation (5),

$\hat{\lambda}_{g}(t|\mathcal{H}_{t})=\sum_{i=1}^{K}\mu(A_{i})\hat{\lambda}_{0i}=\sum_{i=1}^{K}\frac{\mu(A_{i})N_{i}(T)}{T\mu(A_{i})}=\frac{N(T)}{T}$

which is exactly the MLE for the rate of an unmarked Poisson process on the
set $[0,T]$. Hence, $\hat{\lambda}_{g}(t|\mathcal{H}_{t})$ converges in
probability to $\Lambda$. Consider the estimated density function computed
using Equation (6)

$\displaystyle\hat{f}(M)=\frac{\sum_{i=1}^{K}\frac{N_{i}(T)\chi_{A_{i}}(M)}{T\mu(A_{i})}}{\sum_{i=1}^{K}\frac{N_{i}(T)}{T}}=\frac{T}{N(T)}\sum_{i=1}^{K}\frac{N_{i}(T)\chi_{A_{i}}(M)}{T\mu(A_{i})}$
$\implies\hat{f}(M)=\frac{1}{N(T)}\sum_{i=1}^{K}\frac{N_{i}(T)\chi_{A_{i}}(M)}{\mu(A_{i})}.$
(8)

It is apparent that Equation (8) is the histogram estimator for the mark
density of an i.i.d. sample of size $N(T)$. Hence, if the cells $A_{i}$ are
selected correctly this estimator of $f$ is consistent.

######  Remark.

The assumption of setting $\alpha_{ij}=0$ for every $i,j$ may seem unjustified
theoretically. However, this assumption could be justified in practice, by
using a one-sided Kolmogorov-Smirnov test for the uniformity of
$\\{t_{i}\\}_{1\leq i\leq N(T)}$ on $[0,T]$ which may provide insufficient
evidence to suggest temporal clustering of events.

As a result of Example 1, the analytic properties of a histogram estimator are
automatically inherited for our multivariate representation. When
$\mathcal{M}=[0,1]$ we present the main results of [51] Section 6.2 which
suggests how to optimally select the sets $A_{i}$ so that the histogram
estimator converges to the true density. First, the optimal histogram
estimator $\hat{f}$ to the true density $f$ minimises the mean integrated
square error, $\text{MISE}(\hat{f},f)$.

In this case we should select $A_{i}$ as half-open intervals of the form
$A_{i}=\left[\frac{i-1}{K},\frac{i}{K}\right),$ each of uniform width $h=1/K$;
clearly these sets satisfy properties P1-P4. By Theorem 6.11 of [51], under
mild assumptions on $f$, the width of $A_{i}$ that minimises the MISE is
$\mathcal{O}\left(N(T)^{-1/3}\right)$ in which case
$\text{MISE}\left(\hat{f},f\right)=\mathcal{O}\left(N(T)^{-2/3}\right)$.
Therefore, for optimal statistical efficiency in the case of real valued
marks, one should select $K=\lceil N(T)^{-1/3}\rceil$ and uniform bins $A_{i}$
of width $K^{-1}$. For more information see Section 3 of [19] and references
therein.

While analytic results for the MLE given the multivariate representation of a
compound Poisson processes can be derived and shown to be equivalent to known
estimators, history dependent processes are of primary interest. The following
discussion formalises how a multivariate Hawkes representation can approximate
a broad class of non-separable marked Hawkes processes. Approximation of
separable marked Hawkes processes is shown as a special case, as such
processes are subsumed within the larger class of non-separable marked Hawkes
processes.

###  3.2 Non-separable Point Processes

For the following discussion we use a fixed realisation for each process, so
that the conditional intensity functions have discontinuities at the same
instants in time. Underlying our discussion is the result that any measurable
function can be approximated by a simple function, which is the intuition for
Theorems 3.2 and 3.3.

This section is devoted to showing that our multivariate representation is
able to represent a large class of non-separable univariate marked Hawkes
processes. Suppose that $(\mathcal{I}_{1},\mathcal{F}_{1},\mathbb{P}_{1})$ and
$(\mathcal{I}_{2},\mathcal{F}_{2},\mathbb{P}_{2})$ are probability spaces and
that

$f_{1}:\mathcal{I}_{1}\to L^{1}(\mathcal{M}),\quad f_{2}:\mathcal{I}_{2}\to
L^{1}(\mathcal{M}),\quad\xi:\mathcal{I}_{2}\to C(\mathcal{M}).$

Specifically, $f_{1},f_{2}$ are density valued random variables and $\xi$ is a
continuous function valued random variable for all of the mark space
$\mathcal{M}$. Specifically, for every
$\omega_{1}\in\mathcal{I}_{1},\omega_{2}\in\mathcal{I}_{2}$,
$f_{1}(M;\omega_{1})$ and $f_{2}(M;\omega_{2})$ are densities on $\mathcal{M}$
with respect to the same reference measure $\mu$ and $\xi(M;\omega_{2})$ is a
continuous real-valued function known as the mark variable productivity, a
continuous analog to those discussed in [11, 44]. We next consider the general
form of a conditional intensity characteristic of a marked linear Hawkes
process

$\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})=\Lambda\int_{\mathcal{I}_{1}}f_{1}(M;\omega)d\mathbb{P}_{1}(\omega)+\\\
\int_{\mathcal{I}_{2}}f_{2}(M;\omega)\sum_{i:t_{i}<t}g(t-t_{i};M_{i})\xi(M_{i};\omega)d\mathbb{P}_{2}(\omega)$
(9)

where $g(t;M)$ is the time varying kernel dependent on the parameter
$\beta(M)$, such that for every $M\in\mathcal{M}$,
$\int_{0}^{\infty}g(s;M)ds=1$.

The incremental process of arrivals for a given realisation of a non-separable
Hawkes process can be thought of in stages. Since $f_{1},f_{2}$ are densities
pointwise on $\mathcal{I}_{1},\mathcal{I}_{2}$ respectively, their expected
values are as well ($f_{1},f_{2}$ are non-negative so use Tonelli’s theorem to
integrate over $\mathcal{M}$ first). We now apply Equation (5) and find that
points arrive with intensity

$\lambda_{g}(t|\mathcal{H}_{t})=\Lambda+\sum_{i:t_{i}<t}g(t-t_{i})\int_{\mathcal{I}_{2}}\xi(M_{i};\omega)d\mathbb{P}(\omega),$

which is exactly the ground process of a marked Hawkes process. Once a point
arrives, at time $t$, it is an immigrant (in which case set $\tau=1$) with
probability

$p=\frac{\Lambda}{\Lambda+\sum_{i:t_{i}<t}g(t-t_{i})\int_{\mathcal{I}_{2}}\xi(M_{i};\omega)d\mathbb{P}_{2}(\omega)},$

or an offspring event (in which case $\tau=2$) with probability $1-p$. Then
the mark of that point is simulated as a realisation of the random variable
with density $f_{\tau}(M;\omega_{\tau})$. Heuristically, for a given point in
a given realisation, the mark value depends on three rolls of dependent dice.

We now formalise the result, in Theorem 3.2, that there exists a parameter
vector for Equation (4) such that $\lambda_{\theta}(t,M|\mathcal{H}_{t})$ can
approximate $\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})$ well in the $L^{1}$
sense.

######  Theorem 3.2 (Existence).

Consider a non-separable marked Hawkes process with intensity given by
Equation (9) where $f_{1},f_{2}$ are densities with respect to the measure
$\mu$ and $\xi$ is a real-valued continuous function. Moreover, suppose that
$f_{1},f_{2},\xi$ are all equicontinuous on $\mathcal{M}$ with respect to
$\mathcal{I}_{1}$ and $\mathcal{I}_{2}$. Furthermore, assume that $g(t;\beta)$
is uniformly continuous with respect to $\beta$ and that $\beta(M)$ is
continuous in $\mathcal{M}$. Then for every $T>0$, fixed realisation, and
$\epsilon>0$ there exists a $K\in\mathbb{N}$, a parameter vector
$\tilde{\theta}_{K}$, and set of sets $\left\\{A_{i}\right\\}_{1\leq i\leq K}$
that satisfies properties P1-P4 such that,

$\left\lVert\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})-\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})\right\rVert_{L^{1}([0,T]\times\mathcal{M})}<\epsilon.$

###### Proof Sketch.

The sketch of the proof is as follows. We first show

$\bar{f}_{1}(M)=\int_{\mathcal{I}_{1}}f_{1}(M;\omega_{1})d\mathbb{P}_{1}(\omega_{1}),\quad\overline{f_{2}\xi}(M,M^{\prime})=\int_{\mathcal{I}_{2}}f_{2}(M;\omega)\xi(M^{\prime};\omega)\mathbb{P}_{2}(\omega),$

are uniformly continuous on $\mathcal{M}$. We then use compactness to
decompose $\mathcal{M}$ into sets $A_{i}$ so that we can approximate
$\bar{f}_{1}$ and $\overline{f_{2}\xi}$ by a simple function arbitrarily well.
$L^{1}$ convergence is a directly follows. A complete proof is presented in
Appendix A. ∎

We note that Theorem 3.2 could be extended component-wise to a marked
multivariate Hawkes process with $J$ component intensities. Furthermore, in
specific cases where $\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})$ varies
deterministically in time it is possible to select a different parameter
vector $\tilde{\theta}_{K}$ to still approximate
$\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})$, if the representation is selected
wisely. Example 2 demonstrates this for a mark-separable renewal process.

######  Example 2.

Consider a renewal mark-separable Hawkes process, where the immigrant events
arrive as a Weibull renewal process. Hence

$\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})=f(M)\left(ba^{b}(t-t^{\prime})^{b-1}+\sum_{i:t_{i}<t}g(t-t_{i})\xi(M_{i})\right)$

where $a,b$ are parameters and $t^{\prime}$ is the most recent immigrant
event. Then to approximate $\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})$ by our
multivariate representation the parameter vector is then of the form
$\theta=\\{\\{a_{i}\\},\\{b_{i}\\}\\{\alpha_{ij}\\},\\{\beta_{ij}\\}\\}_{{1\leq
i,j\leq K}}$. Pick $\left\\{A_{i}\right\\}_{1\leq i\leq K}$, $\alpha_{ij}$ and
$\beta_{ij}$ as in the proof of Theorem 3.2. Additionally, assume that an
immigrant event resets the renewal rate independent of its mark and set
$b_{i}=b$ and $a_{i}=a\left(\tilde{f}_{i}\right)^{\frac{1}{b_{i}}}$. In this
case, using $\lambda_{0i}(t)=b_{i}a_{i}^{b_{i}}(t-t^{\prime})^{b_{i}-1}$ will
approximate $f(M)ba^{b}(t-t^{\prime})^{b-1}$ in the $L^{1}$ sense.

We now show that the multivariate representation commensurate with the
parameter vector such that $\lambda_{\theta}(t,M)$ approximates a non-
separable Hawkes process on the domain $[0,T]\times\mathcal{M}$ is stationary
if the target Hawkes process is stationary, and for each event its mark-
variable productivity is deterministic.

######  Theorem 3.3 (Stationarity).

Suppose the assumptions of Theorem 3.2 are satisfied and the Hawkes process
$\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})$ is stationary when $\xi$ is
independent of $\omega$. Then
$\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})$, constructed in Theorem
3.2, is also stationary.

###### Proof Sketch.

The stationarity condition for the non-separable Hawkes processes is that the
expected cluster size is finite. To prove this, we first demonstrate that the
expected number of offspring from an immigrant is finite. We then show that
the total number of offspring from a non-immigrant is also finite. In this
case, each independent branching process with one immigrant in the centre,
a.s. has finitely many events. The condition for this is that

$I=\int_{\mathcal{M}}\int_{\mathcal{I}_{2}}f_{2}(M;\omega)\xi(M)d\mathbb{P}_{2}(\omega)d\mu(M)<1.$

Using the parameter vector $\tilde{\theta}_{K}$ and sets $\\{A_{i}\\}_{1\leq
i\leq K}$, as constructed in the proof of Theorem 3.2, we approximate $I$ by
the integral of a simple function. We then compute the expected number of
first generation offspring for our multivariate representation and show that
it is strictly bounded above by $1$ demonstrating that the mean cluster size
is finite a.s. concluding our argument. The full proof is presented in
Appendix B. ∎

######  Remark.

We conjecture Theorem 3.3 to be true for a stochastic mark variable
productivity. This remains conjectured as it is challenging to compute the
expected number of first generation offspring from a multivariate
representation. To do so, Lemma B.1 must be modified so that the mark
distribution is fully history dependent. This in and of itself is not an
issue, but in this case $\alpha_{ij}$ cannot be factorised into a term
dependent on $i$ and a term dependent on $j$. The computation of the expected
number of first generation offspring from a non-immigrant event relies on the
factorisation of $\alpha_{ij}$, and therefore in this context is intractable.
However, in most, if not all, applications $\xi$ is assumed to be
deterministic so this slight loss of generality is not a major concern.

While an expected result, Theorem 3.3 demonstrates that our multivariate
representation has the potential to be used in practice. Of primary
importance, Theorem 3.3 implies the existence of the multivariate
representation $\lambda_{\theta}$ via Khinchin’s Existence Theorem
(Proposition 3.3.I of [12]), as well as orderliness of $\lambda_{\theta}$, a
result upon which the proof of Theorem 3.5 is implicitly reliant [12]. Also of
the utmost importance is that the preservation of stationarity, and hence non-
explosivity almost surely (a.s.), is crucial for frequentist forecasting.
Predictive performance is often performed by repeatedly simulating
realisations of the process over some subset of the time domain [48, 52].
Therefore, if the data generating process is finite, then once the
multivariate representation approximates the data generating process well
enough it too is stationary facilitating the construction of frequentist
forecasts. Finally, the results of Theorem 3.2 imply that our representation
can approximate the true data generating process arbitrarily well.

We again note that mark-separable Hawkes processes are a special case of
Equation (9), which we now discuss in the following example.

######  Example 3.

Consider the mark separable Hawkes process on $[0,T]\times[M_{0},M_{m}]$
induced by the conditional intensity function

$\lambda^{*}(t,M|\mathcal{H}_{t})\vcentcolon=f(M)\left(\Lambda+\sum_{i:t_{i}<t}g(t-t_{i})\xi(M_{i})\right)$
(10)

where $f$ is a continuous density. It is a direct result of Theorems 3.2 and
3.3 that our multivariate representation can approximate $\lambda^{*}$
arbitrarily well, and is stationary if the process induced by $\lambda^{*}$ is
stationary. In the case of the mark separable Hawkes process this is that

$I=\int_{\mathcal{M}}f(M)\xi(M)d\mu(M)<1.$

We now present a corollary, which says that the induced mark density of the
multivariate representation is a simple function approximation of the mark
density $f$.

######  Corollary 3.4.

Suppose the target process is a mark separable Hawkes process with continuous
mark density and mark variable productivity given by Equation (10). Then the
mark density induced by $\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})$ is

$\tilde{f}(M)=\sum_{i=1}^{K}\bar{f}_{i}\chi_{A_{i}}(M),\quad\text{ where
}\bar{f}_{i}=\frac{1}{\mu(A_{i})}\int_{A_{i}}f(M)d\mu(M).$

###### Proof.

Since the conditions of Theorem 3.2 are satisfied we compute the induced mark
density, $f^{\prime}$, using Equation (6) in Proposition 3.1 with
$\tilde{\theta}_{K}$ and the set of sets $\\{A_{i}\\}$ for $1\leq i\leq K$
that exist as a consequence of the theorem. Specifically,
$\lambda_{0i}=\Lambda\bar{f}_{i}$, $\beta_{ij}=\beta$ and
$\alpha_{ij}=\bar{f}_{i}\xi_{j}$. Then

$\displaystyle\tilde{f}(M)=$
$\displaystyle\frac{\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})}{\int_{\mathcal{M}}\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})d\mu(M)}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{K}\chi_{A_{i}}(M)\left(\lambda_{0i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})\right)}{\sum_{i=1}^{K}\mu(A_{i})\left(\lambda_{0i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})\right)}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{K}\chi_{A_{i}}(M)\left(\Lambda\bar{f}_{i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}\bar{f}_{i}g(t-t_{\ell,j})\right)}{\sum_{i=1}^{K}\mu(A_{i})\left(\Lambda\bar{f}_{i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}\bar{f}_{i}g(t-t_{\ell,j})\right)}$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{K}\bar{f}_{i}\chi_{A_{i}}(M)\left(\Lambda+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}g(t-t_{\ell,j})\right)}{\sum_{i=1}^{K}\bar{f}_{i}\mu(A_{i})\left(\Lambda+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}g(t-t_{\ell,j})\right)}=\sum_{i=1}^{K}\bar{f}_{i}\chi_{A_{i}}(M)$

where we obtained the finally equality since
$\sum_{i=1}^{K}\bar{f}_{i}\mu(A_{i})=1$. In particular the marks of the
multivariate representation are i.i.d. for every $K\in\mathbb{N}$, which is as
expected considering the target process has i.i.d. marks. ∎

Corollary 3.4 implies that when modelling potentially separable processes, the
magnitude density can be estimated in a non-parametric way assuming
$\hat{\theta}_{T}\xrightarrow[N(T)\to\infty]{p}\tilde{\theta}_{K}$ as a
consequence of the continuous mapping theorem [25] since
$\lambda_{\theta}(t,M|\mathcal{H}_{t})$ is continuous with respect to
$\theta$. Therefore, once asymptotic distributions for the maximum likelihood
estimates of multivariate processes are established analytically, it may be
possible to develop tests for mark separability and compare it to existing
methods such as [43].

###  3.3 Identifiability

We conclude this section by demonstrating that the parameters of the
multivariate representation are identifiable. Previous studies have
partitioned the mark space in a similar way to us, however they have all
performed non-parametric estimation e.g. [4, 26]. Our multivariate
representation is parametric, and because information is lost by turning
continuous marks (for example if $\mathcal{M}\subset\mathbb{R}^{d}$) into
discrete types, parameter identifiability is not necessarily trivial.
Fortunately, given the below assumptions, the parameters of our multivariate
representation defined by Equation (4), are identifiable.

  * A1.

The point process with intensity $\lambda_{\theta}(t,M|\mathcal{H}_{t})$ is
stationary.

  * A2.

The intensity $\lambda_{\theta}(t,M|\mathcal{H}_{t})$ is positive for every
$(t,M)\in[0,T]\times\mathcal{M}$.

  * A3.

$g_{ij}(t)$ is linearly independent of $g^{\prime}_{ij}(t)$ when
$\beta_{ij}\neq\beta^{\prime}_{ij}$.

  * A4.

In the time interval $[0,T]$ there exists an event such that $(t_{k},M_{k})$
such that $M_{k}\in A_{j}$ for every $j$ and that $t_{1}>0$.

Under assumption A1 the process is finite since it is both simple and
stationary e.g. Proposition 3.3.V of [12]. In this case the process is
orderly, and $\exists\epsilon_{i}>0$, $\forall i\in\mathbb{N}$, such that
$t_{i+1}-t_{i}>\epsilon_{i}$ a.s..

Assumption A2 can be relaxed slightly so that the intensity can be 0 for
strict subsets of $[0,T]\times\mathcal{M}$. However in doing so, we would be
required to define restart times, similar to [6], ensuring no event of type
$i$ arrives at time $t$ when $\lambda_{\theta}(t,M|\mathcal{H}_{t})=0$ for
every $M\in A_{i}$. We leave the close extension of identifiability for
parametric estimation of nonlinear multivariate Hawkes processes
characterising univariate marked processes for future work.

The purpose of assumption A3 is to ensure that for a unique set of parameter
values the conditional intensity function changes in time uniquely. There are
several types of parametric functions that satisfy assumption A3 and have been
used in univariate Hawkes process models. Examples include exponential
functions $g_{ij}(t)=\exp\\{-\beta_{ij}t\\}$ [18, 23], power law functions for
a fixed $p>0$, $g_{ij}(t)=(t+c_{ij})^{-p}$ [31, 32], or
$g_{ij}(t)=te^{-\beta_{ij}t}$ [33, 49].

Finally, assumption A4 ensures that there exists a non-empty subset of
$[0,T]\times\mathcal{M}$ where the intensity depends on every parameter.
Heuristically, if there was never an event of type $1$, i.e. $M_{k}\not\in
A_{1}$ for every $k$, then $\lambda_{\theta}(t,M|\mathcal{H}_{t})$ will never
depend on $\alpha_{11}$. In which case the probability structure of the
process is independent of $\alpha_{11}$, and hence this parameter should not
be identifiable.

######  Theorem 3.5 (Identifiability).

Fix $K\in\mathbb{N}$ to be arbitrary. Assume that the set of sets
$\left\\{A_{j}\right\\}_{1\leq j\leq K}$ satisfy properties P1-P4 and that
assumptions A1-A4 are satisfied. Then parameters of the multivariate
representation, specified by Equation (4), are identifiable.

###### Proof Sketch.

Because conditional intensity functions uniquely define the finite dimensional
distributions of the point process [12], to show identifiability we must
demonstrate that for almost every $(t,M)\in[0,T]\times\mathcal{M}$ that
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\iff\theta=\theta^{\prime}$.
It is obvious that if $\theta=\theta^{\prime}$ then
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})$.

Our argument for the other direction relies on partitioning
$[0,T]\times\mathcal{M}$ into cells of the form $(t_{k},t_{k+1}]\times A_{i}$
(e.g. Figure 1) where $t_{k}$ is the $k^{th}$ event. We then prove
identifiability of the model in the cells $(t_{1},t_{2}]\times A_{i}$, where
$t_{1}$ is the first event. We then induct over the event types to prove the
model is identifiable in the next cell where the conditional intensity
function depends on new parameters. The full proof of Theorem 3.5 is presented
in Appendix C.

0$M_{0}$$t_{1}$$t_{2}$$t_{3}$$A_{4}$$A_{3}$$A_{2}$$A_{1}$$A_{4}$$A_{3}$$A_{2}$$A_{1}$$A_{4}$$A_{3}$$A_{2}$$A_{1}$$A_{4}$$A_{3}$$A_{2}$$A_{1}$
Figure 1: Partitioning of $[0,T]\times\mathcal{M}$

∎

##  4 Simulation study

The following study provides numerical validation for Theorem 3.2. We observe
convergence of simulated realisations of the univariate process observed on
increasingly large time windows to unmarked multivariate processes with an
increasing number of component intensities. We conclude this study with
discussion of the increased model complexity induced by the (possibly) larger
number of parameters necessitated by a multivariate representation relative to
a univariate representation.

###  4.1 Simulation procedure

We consider a stationary univariate marked Hawkes process with a constant
background rate and exponential kernel function, i.e.

$\lambda_{\text{HP}}(t|\mathcal{H}_{t})=\lambda_{0}+\sum_{i:t_{i}<t}\alpha\exp(-\beta(t-t_{i}))$
(11)

which is a special instance of Equation (9) where $\xi(M)=\alpha$ and
$f_{1}=f_{2}$ where both are constant in $\omega$. In this case, the time
varying kernel of the target Hawkes process, $\alpha\exp(-\beta t)$, is not a
density but it could easily be normalised to. Such a process is canonical in
the Hawkes literature, and can be represented as an unmarked univariate Hawkes
process. It was chosen here because the separable uniform mark distribution is
maximally entropic, and therefore minimally informative with respect to the
ground intensity. This means that the estimates of the multivariate Hawkes
representation will be maximally statistically inefficient, and therefore this
study provides a heuristic “upper bound” on the error induced by the increased
number of parameters necessitated by the multivariate representation. Marks
were simulated as draws from a discrete uniform distribution
$f\sim\text{Unif}(1,K)$.

We simulated S=128 realisations on the time window $[0,5056]$. Subsets of each
realisation were then observed on twenty increasingly large time windows
commensurate with a range of $10^{2}$ to $10^{4}$ points. Maximum likelihood
parameter estimates were obtained for multivariate Hawkes process
representations $\lambda_{\theta}(t,M|\mathcal{H}_{t})$ given $K=1,2,\ldots,6$
component intensities and an exponential kernel matrix. The properties of such
estimates are discussed in [36]; consistency of MLE for multivariate Hawkes
processes is widely conjectured and numerically validated [7]. The parameter
vector $\theta^{*}=\\{\alpha^{*}=1,\beta=^{*}2,\mu^{*}=1\\}$ was specified as
$\int_{0}^{\infty}\alpha^{*}\exp(-\beta^{*}u)du=\frac{1}{2}$, guaranteeing
stationarity.

###  4.2 Simulation results and interpretation

Given specified parameter values and a uniform mark distribution, ansatz
parameter values are equal to
$\bm{\theta}^{*}=\\{\bm{\alpha}^{*},\bm{\beta}^{*},\bm{\mu}^{*}\\}$ where the
vector $\bm{\mu}=\mu^{*}K^{-1}\bm{1}_{K}$,
$\bm{\alpha}=\alpha^{*}K^{-1}\bm{1}_{K}\bm{1}_{K}^{\top}$ and
$\bm{\beta}={\beta^{*}}\bm{1}_{K}\bm{1}_{K}^{\top}$ where
$\bm{1}_{K}=(1,1,\ldots,1)\in\mathbb{R}^{K}$. Median absolute error (MAE) was
then computed for each MLE estimated parameter vector $\bm{\hat{\theta}}$,
defined by
$\text{MAE}(\\{\bm{\hat{\theta}}_{s}\\}_{s=1}^{S},\bm{\theta}^{*})=\text{median}\\{||\bm{\hat{\theta}}_{s}-\bm{\theta}^{*}||_{1}\\}_{s=1}^{S}$.

Figure 2: Error associated with parameter estimates of multivariate Hawkes
representations of univariate marked Hawkes process with an exponential kernel
function. Plots show median $L^{1}$ error of parameter estimates across 128
realisations of size ranging from $10^{2}$ to $10^{4}$ points. Error is
generally larger as the number of components $K$ increases. $90$% confidence
intervals are shown in appendix D. Number of points and error are shown on the
$\log_{10}$ scale.

Multivariate unmarked Hawkes process representations of univariate marked
Hawkes processes have more parameters in most cases. For example, given an
exponential kernel matrix and constant background parameters, a multivariate
Hawkes process has with $K$ components, $K+2K^{2}$ parameters. A marked
univariate Hawkes process with exponential kernel has $3+M$ parameters where
$M$ is the dimension of the parameter space of the mark distribution, and in
most cases, $K+2K^{2}>3+M$. However, despite the larger number of parameters,
given its flexibility to characterise mark distributions quasi-parametrically,
the multivariate representation is a arguably the more parsimonious and
computationally tractable model in many scenarios of interest.

A discrete uniform mark distribution was specified because it is maximally
entropic (i.e. minimally informative with respect to the ground intensity) and
therefore, the results visualized in Figure 2 can be used as a rough heuristic
bound for worst-case MAE induced by the larger number of parameters necessary
to specify a multivariate Hawkes representation of a univariate marked
process. In particular, simulated MAE appears $\mathcal{O}\left(K\cdot
N(T)^{-\frac{1}{2}}\right)$.

Comparison with the cLSTM neural network architecture [27] with regards to the
number of parameters required by the multivariate representation is most apt
as the cLSTM allows for univariate integer-marked Hawkes process data to be
represented by an unmarked parametric multivariate Neural Hawkes process. The
cLSTM architecture has the fewest parameters relative to other options such as
a fully connected feed forward architecture [35] or transformer architecture
[55], but still has $d(2K+14d+1)$ parameters given $d$ hidden units and $K$
component intensities. The need for regularization of multivariate Hawkes
processes is well discussed in the literature, for example see [2], and we
leave the ramifications of sparsity inducing regularization in the context of
multivariate Hawkes representations of univariate marked Hawkes processes to
future work. All code for this study is freely available at
https://github.com/ckres213/MVH_Rep.

##  5 Conclusion

In this article we propose a multivariate Hawkes representation for
parameterizing univariate marked Hawkes processes. We demonstrate its
stationarity, identifiability, and that for a given realisation, the
conditional intensity functions of the representative and approximated
processes converge in $L^{1}$ as the number of components increases. We
further present a simulation study demonstrating the $L^{1}$ convergence of
parameter estimates to theoretical values given a maximally entropic mark
distribution.

We suggest two possible avenues for future research. Because the density
estimator for our process is akin to a histogram estimator, which is clear
from Example 1, we hypothesise a setup which is analogous to that of a kernel
density estimator (see [5] for an application to spatial point processes).
This would likely increase the ease of application for the multivariate
representations, especially in the case of high dimensional mark spaces.
Finally, inducing sparsity to reduce the dimension of the parameter space by
coalescing subcells of the mark space would likely decrease the model
complexity of the representation, thereby increasing its usability.

##  Appendix A Proof of Theorem 3.2

###### Proof.

We show that for each $T>0$ and $\epsilon>0$ there exists a $K\in\mathbb{N}$,
parameter vector $\tilde{\theta}_{K}$, and sets $\left\\{A_{i}\right\\}_{1\leq
i\leq K}$ satisfying properties P1-P4, such that

$I\vcentcolon=\left\lVert\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})-\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})\right\rVert_{L^{1}([0,T]\times\mathcal{M})}<\epsilon.$
(12)

We first present a relevant lemma.

######  Lemma A.1.

Suppose that $(\mathcal{I}_{1},\mathcal{F}_{1},\mathbb{P}_{1})$ and
$(\mathcal{I}_{2},\mathcal{F}_{2},\mathbb{P}_{2})$ are probability spaces and
that

$f_{1}:\mathcal{I}_{1}\to L^{1}(\mathcal{M}),\quad f_{2}:\mathcal{I}_{2}\to
L^{1}(\mathcal{M}),\quad\xi:\mathcal{I}_{2}\to C(\mathcal{M})$

where $f_{1},$ $f_{2}$ and $\xi$ are equicontinuous on $\mathcal{M}$, a
compact CSMS, for $\mathcal{I}_{1}$ or $\mathcal{I}_{2}$ respectively. Define

$\bar{f}_{1}(M)=\int_{\mathcal{I}_{1}}f_{1}(M;\omega_{1})d\mathbb{P}_{1}(\omega_{1}),\quad\overline{f_{2}\xi}(M,M^{\prime})=\int_{\mathcal{I}_{2}}f_{2}(M;\omega)\xi(M^{\prime};\omega)d\mathbb{P}_{2}(\omega),$

then $\bar{f}_{1}$ and $\overline{f_{2}\xi}(M,M^{\prime})$ are uniformly
continuous on $\mathcal{M}$ and $\mathcal{M}^{2}$ respectively.

###### Proof.

Since $\mathcal{M}$ is compact $f_{1}(\cdot;\omega)$ is uniformly continuous
for each $\omega\in\mathcal{I}_{i}$. Moreover, $\mathcal{M}^{2}$ is compact so
$f_{2}(M;\omega)\xi(M^{\prime};\omega)$ is also uniformly continuous for each
$\omega\in\mathcal{I}_{2}$. Hence for every $\epsilon>0$,
$\omega\in\mathcal{I}_{1}$ and every $M\in\mathcal{M}$, there exists a
$\delta>0$ so that if $d(M,M^{\prime})<\delta$ then
$|f_{1}(M,\omega)-f_{1}(M^{\prime},\omega)|<\epsilon$. Clearly

$\displaystyle\left|\bar{f}_{1}(M)-\bar{f}_{1}(M^{\prime})\right|\leq\int_{\mathcal{I}_{1}}\left|f_{1}(M;\omega)-f_{1}(M^{\prime};\omega)\right|d\mathbb{P}_{1}(\omega)<\epsilon$

since $\mathbb{P}$ is a probability measure. This exact argument can be
extended to $\overline{f_{2}\xi}(M,M^{\prime})$. ∎

We first construct the sets $\left\\{A_{i}\right\\}_{1\leq i\leq K}$ that
cover $\mathcal{M}$ and satisfy properties P1-P4. For any $\epsilon>0$, cover
$\mathcal{M}$ by the finitely many balls $\mathbb{B}(M_{i},\delta)$ such that
for any $M,M^{\prime}\in\mathbb{B}(M_{i},\delta)$ and
$m,m^{\prime}\in\mathbb{B}(M_{j},\delta)$ we have

$|\bar{f}_{1}(M)-\bar{f}_{1}(M^{\prime})|<\epsilon,\quad|\overline{f_{2}\xi}(M,m)-\overline{f_{2}\xi}(M^{\prime},m^{\prime})|<\epsilon,\quad|g(t,M)-g(t,M^{\prime})|<\epsilon\quad\forall
t.$

We now construct

$A_{i}=\left(\mathbb{B}(M_{i},\delta)\setminus\bigcap_{j=1}^{i-1}\mathbb{B}(M_{j},\delta)\right)$

and discard the sets that have measure $0$. Since we assume $\mathcal{M}$ has
finite measure, then so too does every $A_{i}\subset\mathcal{M}$. Clearly,
properties P1-P4 are satisfied by these sets.

Now we unpack definitions and notation to find

$I=\int_{0}^{T}\int_{\mathcal{M}}\Bigg{|}\bar{f}_{1}(M)\Lambda+\sum_{\ell:t_{\ell}<t}g(t-t_{\ell};M_{\ell})\overline{f_{2}\xi}(M,M_{\ell})-\\\
\sum_{i=1}^{K}\chi_{A_{i}}(M)\left(\lambda_{0i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})\right)\Bigg{|}d\mu(M)dt.$

By construction $\bigcup_{i=1}^{K}A_{i}$ covers $\mathcal{M}$ and $A_{m}\cap
A_{n}=\emptyset$ when $m\neq n$ therefore

$I=\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\Bigg{|}\bar{f}_{1}(M)\Lambda+\sum_{\ell:t_{\ell}<t}g(t-t_{\ell};M_{\ell})\overline{f_{2}\xi}(M,M_{\ell})-\\\
\left(\lambda_{0i}+\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})\right)\Bigg{|}d\mu(M)dt$

$\leq\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\left|\bar{f}_{1}(M)\Lambda-\lambda_{0i}\right|d\mu(M)dt+\\\
\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\left|\sum_{\ell:t_{\ell}<t}g(t-t_{\ell};M_{\ell})\overline{f_{2}\xi}(M,M_{\ell})-\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})\right|d\mu(M)dt.$

Denote

$\tilde{f}_{i,1}=\frac{1}{\mu(A_{i})}\int_{A_{i}}\bar{f}_{1}(M)d\mu(M),\quad\widetilde{f_{2}\xi}_{ij}=\frac{1}{\mu(A_{i}\times
A_{j})}\int_{A_{i}\times A_{j}}\overline{f_{2}\xi}(M,M^{\prime})d\mu(M\times
M^{\prime}).$

Define $\lambda_{0i}=\Lambda\tilde{f}_{i,1}$,
$\beta_{ij}\in\beta(A_{j})=:\\{\beta(M)\in\mathbb{R}:M\in A_{j}\\}$, and pick
$\alpha_{ij}=\widetilde{f_{2}\xi}_{ij}$. Re-index the middle sum to
$\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}$ which is summing over each
event once by type, as opposed to their arrival time. Then

$I\leq\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\Lambda\left|\bar{f}_{1}(M)-\tilde{f}_{i,1}\right|d\mu(M)dt+\\\
\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\left|\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}g(t-t_{\ell,j};M_{\ell,j})\overline{f_{2}\xi}(M,M_{\ell,j})-\widetilde{f_{2}\xi}_{ij}g_{ij}(t-t_{\ell,j})\right|d\mu(M)dt.$

$\implies
I\leq\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\Lambda\left|\bar{f}_{1}(M)-\tilde{f}_{i,1}\right|d\mu(M)dt+\\\
\hskip
48.36958pt\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\overline{f_{2}\xi}(M,M_{\ell,j})\left|g(t-t_{\ell,j};M_{\ell,j})-g_{j}(t-t_{\ell,j})\right|d\mu(M)dt+\\\
\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}g_{ij}(t-t_{\ell,j})\left|\overline{f_{2}\xi}(M,M_{\ell,j})-\widetilde{f_{2}\xi}_{ij}\right|d\mu(M)dt\\\
=\vcentcolon\text{LHS}+\text{MT}+\text{RHS}.$ (13)

Consider the first term in Equation (13) (LHS). We change the order of
integration and integrate over $t$ first. Since $\bar{f}_{1}$ is uniformly
continuous $\left|\bar{f}_{1}(M)-\tilde{f}_{i,1}\right|<\epsilon$. Hence,
$\text{LHS}<\mu(\mathcal{M})\Lambda T\epsilon$.

Before we bound the next term (MT), we define
$\bm{f\xi}\vcentcolon=\sup_{(M,M^{\prime})\in\mathcal{M}^{2}}\overline{f_{2}\xi}(M,M_{\ell,j})<\infty$
since $\overline{f_{2}\xi}(M,M_{\ell,j})$ is continuous and $\mathcal{M}^{2}$
is compact. Then using the continuity of $g(t;M)$, with respect to $M$,

$\displaystyle\text{MT}\leq$
$\displaystyle\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\bm{f\xi}\left|g(t-t_{\ell,j};M_{\ell,j})-g_{j}(t-t_{\ell,j})\right|d\mu(M)dt$
$\displaystyle<$
$\displaystyle\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\bm{f\xi}\epsilon
d\mu(M)dt$ $\displaystyle<$ $\displaystyle\bm{f\xi}\epsilon
N(T)\mu(\mathcal{M}).$

To bound the final term, we first use continuity to find that
$\left|\overline{f_{2}\xi}(M,M_{\ell,j})-\widetilde{f_{2}\xi}_{ij}\right|<\epsilon.$
We also use that for every $i,j\in\\{1,2,\dots K\\}$, $g_{ij}(t)$ is a density
with respect to $t$. Therefore,

$\displaystyle RHS\leq$
$\displaystyle\int_{0}^{T}\sum_{i=1}^{K}\int_{A_{i}}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}g_{ij}(t-t_{\ell,j})\epsilon
d\mu(M)dt$ $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{K}\int_{A_{i}}N(T)\epsilon
d\mu(M)=\mu(\mathcal{M})N(T)\epsilon.$

Hence, for every $\varepsilon>0$ we can find a $K\in\mathbb{N}$, sets
$\left\\{A_{i}\right\\}_{1\leq i\leq K}$ and parameter vector
$\tilde{\theta}_{K}$ such that

$\left\lVert\lambda_{\text{HP}}(t,M|\mathcal{H}_{t})-\lambda_{\tilde{\theta}}(t,M|\mathcal{H}_{t})\right\rVert_{L^{1}([0,T]\times\mathcal{M})}<\varepsilon.$

By carefully selecting the parameter vector, we can remove the implicit
dependence of $K$ on $T$. This is so that can $K$ can be selected to be
arbitrarily large for any fixed $T$ so that Equation (12) holds for every
finite positive $T$ (i.e. $I\to 0$ as $K\to\infty$ for each fixed $T>0$).

∎

##  Appendix B Proof of Theorem 3.3

###### Proof.

Because the targeted non-separable Hawkes process is stationary, the expected
number of offspring from a non-immigrant event is less than $1$ [54].
Specifically, $I\vcentcolon=\mathbb{E}[\xi(M)|M\text{ has a parent}]<1$ which
can be expressed as

$I=\int_{\mathcal{M}}\left(\int_{\mathcal{I}_{2}}f_{2}(M;\omega)d\mathbb{P}_{2}(\omega)\right)\xi(M)d\mu(M)=1-\sigma<1$

for $\sigma\in(0,1)$. Since the assumptions of Theorem 3.2 hold we use results
from its proof, in the special case that $\xi$ is independent of $\omega$. In
which case, for any $\epsilon>0$ there exists a $K\in\mathbb{N}$ and set of
sets $\\{A_{i}\\}_{1\leq i\leq K}$ satisfying P1-P4 such that for any
$M,M^{\prime}\in A_{i}$,

$\quad|\bar{f}_{2}(M)-\bar{f}_{2}(M^{\prime})|<\epsilon,\quad|\xi(M)-\xi(M^{\prime})|<\epsilon$

where
$\bar{f}_{2}(M)=\int_{\mathcal{I}_{2}}f_{2}(M;\omega)d\mathbb{P}_{2}(\omega)$,
for every $i\in\\{1,2,\dots K\\}$. Then by picking
$\xi_{i}\in\\{\xi(M)\in\mathbb{R}^{+}:M\in A_{i}\\}$ and
$\tilde{f}_{i,2}=\int_{\mathcal{M}}\bar{f}_{2}(M)d\mu(M)$ we examine the
quantity

$\displaystyle\left|\int_{\mathcal{M}}\bar{f}_{2}(M)\xi(M)d\mu(M)-\sum_{i=1}^{K}\tilde{f}_{i,2}\xi_{i}\mu(A_{i})\right|\leq\sum_{i=1}^{K}\int_{A_{i}}\left|\bar{f}_{2}(M)\xi(M)-\tilde{f}_{i,2}\xi_{i}\right|d\mu(M)$
$\displaystyle\leq\sum_{i=1}^{K}\int_{A_{i}}\bar{f}_{2}(M)\left|\xi(M)-\xi_{i}\right|d\mu(M)+\sum_{i=1}^{K}\int_{A_{i}}\xi_{i}\left|\bar{f}_{2}(M)-\tilde{f}_{i,2}\right|d\mu(M)$
$\displaystyle<\epsilon(1+\bar{\xi}\mu(\mathcal{M}))=\vcentcolon\varepsilon$

where $\bar{\xi}\vcentcolon=\sup_{M\in\mathcal{M}}\xi(M)<\infty$ since
$\mathcal{M}$ is compact and $\xi$ is continuous. Thus, for every
$\varepsilon>0$ there exists a $K\in\mathbb{N}$ such that

$J\vcentcolon=\sum_{i=1}^{K}\tilde{f}_{i,2}\xi_{i}\mu(A_{i})<I+\varepsilon=1-\sigma+\varepsilon.$
(14)

In particular if $\varepsilon<\sigma$ then $J<1$.

It is clear that in the case $\xi$ is independent of $\omega$ that the
parameter vector $\tilde{\theta}_{K}$ has entries of the form
$\lambda_{0i}=\Lambda\tilde{f}_{i,1}$, $\beta_{ij}=\beta_{j}$,
$\alpha_{ij}=\tilde{f}_{i,2}\xi_{j}$. Before, presenting the stationarity
condition for our multivariate representation we demonstrate another lemma.

######  Lemma B.1.

For every $K\in\mathbb{N}$ the mark density of the offspring events induced by
$\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})$ is

$f(M)=\sum_{i=1}^{K}\tilde{f}_{i,2}\chi_{A_{i}}(M).$

###### Proof.

For any Borel set $A\in\mathcal{M}$, by Bayes’ Theorem

$\displaystyle\mathbb{P}[m\in A|m\text{ is an
offspring}]=\frac{\mathbb{P}[m\text{ is an offspring}|m\in A]\mathbb{P}[m\in
A]}{\mathbb{P}[m\text{ is an offspring}]}$ $\displaystyle=$
$\displaystyle\frac{\int_{A}\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})d\mu(M)\frac{\int_{A}\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})d\mu(M)}{\int_{\mathcal{M}}\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})d\mu(M)}}{\int_{A}\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})d\mu(M)}\times$
$\displaystyle\left(\frac{\int_{\mathcal{M}}\lambda_{\tilde{\theta}_{K}}(t,M|\mathcal{H}_{t})d\mu(M)}{\int_{\mathcal{M}}\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})d\mu(M)}\right)$
$\displaystyle=$
$\displaystyle\frac{\int_{A}\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})d\mu(M)}{\int_{\mathcal{M}}\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})d\mu(M)}.$

Clearly, the density $f$ with respect to the measure $\mu$ is then

$f(M)=\frac{\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})}{\int_{\mathcal{M}}\sum_{i,j=1}^{K}\chi_{A_{i}}(M)\sum_{t_{\ell,j}:t_{\ell,j}<t}\alpha_{ij}g_{ij}(t-t_{\ell,j})d\mu(M)}.$

Now using the parameter values $\tilde{\theta}_{K}$

$\displaystyle f(M)=$
$\displaystyle\frac{\sum_{i=1}^{K}\chi_{A_{i}}\tilde{f}_{i,2}(M)\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}g_{j}(t-t_{\ell,j})}{\sum_{i=1}^{K}\mu(A_{i})\tilde{f}_{i,2}\sum_{j=1}^{K}\sum_{t_{\ell,j}:t_{\ell,j}<t}\xi_{j}g_{j}(t-t_{\ell,j})}=\sum_{i=1}^{K}\tilde{f}_{i,2}\chi_{A_{i}}(M)$

since $\sum_{i=1}^{K}\mu(A_{i})\tilde{f}_{i,2}=1$. ∎

We now demonstrate that the expected cluster size is finite, hence showing
stationarity for the multivariate representation.

The expected number of first generation offspring of an event of type $j$,
labelled $E_{j}$, is the integral over $[0,T]\times\mathcal{M}$ of its
contribution to the intensity. Therefore,

$E_{j}=\int_{0}^{T}\int_{\mathcal{M}}\sum_{i=1}^{K}\chi_{A_{i}}(M)\alpha_{ij}g(t)d\mu(M)dt=\sum_{i=1}^{K}\mu(A_{i})\alpha_{ij}.$

Alternatively, the expected number of first generation offspring of an event
of size $M$ is

$E(M)=\sum_{j=1}^{K}\chi_{A_{j}}(M)E_{j}=\sum_{j=1}^{K}\chi_{A_{j}}(M)\sum_{i=1}^{K}\mu(A_{i})\alpha_{ij}.$

The expected number of first generation offspring from a non-immigrant event
is then

$\displaystyle E_{C}=$
$\displaystyle\int_{\mathcal{M}}E(M)f(M)d\mu(M)=\int_{\mathcal{M}}\sum_{j=1}^{K}\chi_{A_{j}}(M)\left(\sum_{i=1}^{K}\mu(A_{i})\alpha_{ij}\right)f(M)d\mu(M)$
$\displaystyle=$
$\displaystyle\int_{\mathcal{M}}\sum_{j=1}^{K}\chi_{A_{j}}(M)\sum_{i=1}^{K}\mu(A_{i})\alpha_{ij}\tilde{f}_{j,2}d\mu(M)=\sum_{j=1}^{K}\mu(A_{j})\sum_{i=1}^{K}\mu(A_{i})\tilde{f}_{i,2}\xi_{j}\tilde{f}_{j,2}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{K}\mu(A_{j})\xi_{j}\tilde{f}_{j,2}\sum_{i=1}^{K}\mu(A_{i})\tilde{f}_{i,2}=\sum_{j=1}^{K}\mu(A_{j})\xi_{j}\tilde{f}_{j,2}=J<1-\sigma+\varepsilon.$

Hence we pick $K$ large enough so that $J<1$, which implies $E_{C}<1$.
Therefore, the mean number of first generation offspring of a non-immigrant
event is less than one. It follows that the geometric series equal to the mean
number of total offspring from a non-immigrant event,
$\sum_{r=0}^{\infty}E_{C}^{r}$, is finite. Immigrant events produce a finite
number of offspring events a.s. because

$\int_{\mathcal{M}}\bar{f}_{1}(M)\xi(M)d\mu(M)<\mu(\mathcal{M})\sup_{M\in\mathcal{M}}\bar{f}_{1}(M)\xi(M)<\infty$

by the continuity of the integrand as well as the compactness and finite
measure of $\mathcal{M}$. Then immigrant events a.s. produce only finitely
many first generation offspring, each of which only produce finitely many
offspring in total. Therefore, a.s. every cluster is off a finite size. Since
branching processes are stationary if their mean cluster size is finite (e.g.
Exercise 6.3.5 of [12]) this concludes the proof.

∎

##  Appendix C Proof of Theorem 3.5

###### Proof.

Since conditional intensities uniquely define the finite dimensional
distributions of the point process [12], to show identifiability we want to
demonstrate that for almost every $(t,M)\in[0,T]\times\mathcal{M}$ that
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\iff\theta=\theta^{\prime}$.
First of all, it is clear that if $\theta=\theta^{\prime}$ then
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})$.

We now assume that
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})$
almost everywhere (a.e.) on $[0,T]\times\mathcal{M}$. First consider
$(t,M)\in[0,t_{1}]\times A_{i}$. Since we assume $t_{1}>0$ we find

$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\implies\lambda_{0i}=\lambda_{0i}^{\prime}$

since the history is empty. $i$ was arbitrary, hence $\lambda_{0i}$ is
identifiable for every $i\in\\{1,2,\dots K\\}$.

Consider the first event $t_{1}$ and, by relabelling to ease notation, suppose
that it is of type $1$, i.e. $M_{1}\in A_{1}$. Then for
$(t,M)\in(t_{1},t_{2}]\times A_{i}$ we have

$\displaystyle\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{0i}+\alpha_{i1}g_{i1}(t-t_{1}),\quad\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})=\lambda_{0i}+\alpha^{\prime}_{i1}g^{\prime}_{i1}(t-t_{1})$
$\displaystyle\lambda_{\theta}(t,M|\mathcal{H}_{t})=$
$\displaystyle\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\implies\alpha_{i1}g_{i1}(t-t_{1})=\alpha^{\prime}_{i1}g^{\prime}_{i1}(t-t_{1})$
(15) $\displaystyle\iff$
$\displaystyle\alpha_{i1}g_{i1}(t-t_{1})-\alpha^{\prime}_{i1}g^{\prime}_{i1}(t-t_{1})=0\mbox{
a.e. in }(t_{1},t_{2}]\times A_{i}.$ (16)

By linear independence Equation (16) can only be true if
$g_{i1}(t-t_{1})=g^{\prime}_{i1}(t-t_{1})\iff\beta_{i1}=\beta^{\prime}_{i1}$
or $\alpha_{i1}=\alpha^{\prime}_{i1}=0$.

In the first case, i.e. if $g_{i1}(t-t_{1})=g^{\prime}_{i1}(t-t_{1})$ then

$\alpha_{i1}g_{i1}(t-t_{1})=\alpha^{\prime}_{i1}g^{\prime}_{i1}(t-t_{1})\iff\alpha_{i1}=\alpha^{\prime}_{i1}$

since $g_{i1}(t-t_{1})>0$ by assumption.

Otherwise, $\alpha_{i1}=\alpha^{\prime}_{i1}=0$. This means that $\beta_{i1}$
is unidentifiable, however the probability structure of the process is
independent of $\beta_{i1}$ so it is trivially unidentifiable. Hence, for
every $i$, when the probability structure of the point process depends on
$\alpha_{i1}$ and $\beta_{i1}$ these parameters are identifiable. Hence our
inductive base case is true.

Now suppose we have observed events of type $1,2,\dots L$, as in there exists
an event $(t_{\ell},M_{\ell})$ such that $M_{\ell}\in A_{n}$ for every
$n\in\\{1,2,\dots L\\}$, which we know to be true by assumption A4.

We inductively assume $\alpha_{in}=\alpha^{\prime}_{in}$ and
$\beta_{in}=\beta^{\prime}_{in}$ for $L<K$ and every $i\in\\{1,2,\dots K\\}$.
We now show that $\alpha_{i(L+1)}=\alpha^{\prime}_{i(L+1)}$ and
$\beta_{i(L+1)}=\beta^{\prime}_{i(L+1)}$.

Suppose that $t_{m}$ is the first event of type $L+1$ and that
$\mathcal{K}_{t}$ is the index set of every type of event that has occurred
prior to time $t$. Consider, $(t,M)\in(t_{m},t_{m+1}]\times A_{i}$ where
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})$
a.e., then

$\lambda_{0i}+\sum_{k\in\mathcal{K}_{t}}\sum_{t_{\ell,k}:t_{\ell,k}<t}\alpha_{ik}g_{ik}(t-t_{\ell,k})=\lambda^{\prime}_{0i}+\sum_{k\in\mathcal{K}_{t}}\sum_{t_{\ell,k}:t_{\ell,k}<t}\alpha^{\prime}_{ik}g^{\prime}_{ik}(t-t_{\ell,k})\mbox{
a.e.}.$ (17)

However, we have shown $\lambda_{0i}=\lambda^{\prime}_{0i}$ and by our
inductive hypothesis for every $k\in\mathcal{K}_{t_{m}}=\\{1,2,\dots L\\}$,
and $i\in\\{1,2,\dots K\\}$, $\alpha_{ik}=\alpha^{\prime}_{ik}$ and
$\beta_{ik}=\beta^{\prime}_{ik}$. By cancelling these common terms from both
sides of Equation (17) we find

$\alpha_{i(L+1)}g_{i(L+1)}(t-t_{m})=\alpha^{\prime}_{i(L+1)}g^{\prime}_{i(L+1)}(t-t_{m})\\\
\iff\alpha_{i(L+1)}g_{i(L+1)}(t-t_{m})-\alpha^{\prime}_{i(L+1)}g^{\prime}_{i(L+1)}(t-t_{m})=0.$
(18)

The right hand side of Equation (18) is identical to Equation (16) except one
must replace $1$ by $L+1$. As in this previous case, we find
$\alpha_{i(L+1)}=\alpha^{\prime}_{i(L+1)}$ and
$\beta_{i(L+1)}=\beta^{\prime}_{i(L+1)}$ when the probability structure of the
process depends on these parameters. Therefore, if the $L^{th}$ case is true
then the $(L+1)^{th}$ case is true, i.e. $\alpha_{in}=\alpha^{\prime}_{in}$
and $\beta_{in}=\beta^{\prime}_{in}$ for
$n\in\mathcal{K}_{t_{M+1}}=\\{1,2,\dots L+1\\}$ and every $i\in\\{1,2,\dots
K\\}$. By the principle of mathematical induction we find that
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\mbox{
a.e. }\implies\theta=\theta^{\prime}$.

Therefore,
$\lambda_{\theta}(t,M|\mathcal{H}_{t})=\lambda_{\theta^{\prime}}(t,M|\mathcal{H}_{t})\mbox{
a.e. }\iff\theta=\theta^{\prime}$. Equivalently, all parameters that the
probability structure of the process depends on are identifiable.

∎

##  Appendix D Figures of Section 4

Figure 3: Approximate 90% confidence interval for median $L^{1}$ error
associated with the parameters of the multivariate Hawkes representation of
the univariate Hawkes process specified in Equation 11. Median and associated
confidence interval based on 128 realisations.

##  Appendix E Acknowledgements and Funding

The authors wish to acknowledge the use of New Zealand eScience Infrastructure
(NeSI) high performance computing facilities, and consulting support as part
of this research. New Zealand’s national facilities are provided by NeSI and
funded jointly by NeSI’s collaborator institutions and through the Ministry of
Business, Innovation & Employment’s Research Infrastructure programme.
https://www.nesi.org.nz.

This project was supported by the Royal Society of New Zealand Marsden Fund
(contact MFP-UOO2323).

## References

  * [1] L. Adamopoulos and A. G. Hawkes.  Cluster models for earthquakes-regional comparisons.  Bull. Int. Stat. Inst., 45(3):454–461, 1973.
  * [2] E. Bacry, M. Bompaire, S. Gaïffas, and J.-F. Muzy.  Sparse and low-rank multivariate Hawkes processes.  Journal of Machine Learning Research, 21(50):1–32, 2020.
  * [3] E. Bacry, I. Mastromatteo, and J.-F. Muzy.  Hawkes processes in finance.  Market Microstructure and Liquidity, 1(01):1550005, 2015.
  * [4] E. Bacry and J.-F. Muzy.  First and second order statistics characterization of Hawkes processes and non-parametric estimation.  IEEE Transactions on Information Theory, 62(4):2184–2202, 2016.
  * [5] A. Baddeley, E. Rubak, and R. Turner.  Spatial point patterns: Methodology and applications with R.  CRC press, 2015.
  * [6] A. Bonnet, M. Martinez Herrera, and M. Sangnier.  Inference of multivariate exponential Hawkes processes with inhibition and application to neuronal activity.  Statistics and Computing, 33(4):91, 2023.
  * [7] C. G. Bowsher.  Modelling security market events in continuous time: intensity based, multivariate point process models.  Journal of Econometrics, 141(2):876–912, 2007.
  * [8] P. Brémaud.  Point processes and queues.  Springer, 1981.
  * [9] P. Brémaud and L. Massoulié.  Stability of nonlinear Hawkes processes.  The Annals of Probability, 24(3):1563–1588, 1996.
  * [10] F. Cheysson and G. Lang.  Spectral estimation of Hawkes processes from count data.  The Annals of Statistics, 50(3):1722–1746, 2022.
  * [11] W.-H. Chiang, X. Liu, and G. Mohler.  Hawkes process modeling of COVID-19 with mobility leading indicators and spatial covariates.  International journal of forecasting, 38(2):505–520, 2022.
  * [12] D. J. Daley and D. Vere-Jones.  An introduction to the theory of point processes: volume I: elementary theory and methods.  Springer, 2003.
  * [13] D. J. Daley and D. Vere-Jones.  An introduction to the theory of point processes: volume II: general theory and structure.  Springer, 2008.
  * [14] L. Davis, B. Baeumer, and T. Wang.  A multidimensional fractional Hawkes process for multiple earthquake mainshock aftershock sequences.  arXiv preprint arXiv:2404.01478, 2024.
  * [15] P. Embrechts and M. Kirchner.  Hawkes graphs.  Theory of Probability & Its Applications, 62(1):132–156, 2018.
  * [16] P. Embrechts, T. Liniger, and L. Lin.  Multivariate Hawkes processes: an application to financial data.  Journal of Applied Probability, 48(A):367–378, 2011.
  * [17] A. S. Fotheringham and D. W. Wong.  The modifiable areal unit problem in multivariate statistical analysis.  Environment and Planning A, 23(7):1025–1044, 1991.
  * [18] A. G. Hawkes.  Spectra of some self-exciting and mutually exciting point processes.  Biometrika, 58(1):83–90, 1971.
  * [19] A. J. Izenman.  Review papers: Recent developments in nonparametric density estimation.  Journal of the American Statistical Association, 86(413):205–224, 1991.
  * [20] J. Jacod.  Multivariate point processes: predictable projection, Radon-Nikodym derivatives, representation of martingales.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 31(3):235–253, 1975.
  * [21] M. Jun and S. Cook.  Flexible multivariate spatiotemporal Hawkes process models of terrorism.  The Annals of Applied Statistics, 18(2):1378 – 1403, 2024.
  * [22] C. Kresin and F. Schoenberg.  Parametric estimation of spatial–temporal point processes using the Stoyan–Grabarnik statistic.  Annals of the Institute of Statistical Mathematics, 75(6):887–909, 2023.
  * [23] C. Kresin, F. P. Schoenberg, and G. Mohler.  Comparison of Hawkes and SEIR models for the spread of covid-19.  Advances and Applications in Statistics, 74:83–106, 2022.
  * [24] T. Liniger.  Multivariate Hawkes processes.  PhD thesis, ETH Zurich, 2009.
  * [25] H. B. Mann and A. Wald.  On stochastic limit and order relationships.  The Annals of Mathematical Statistics, 14(3):217 – 226, 1943.
  * [26] E. B. Marcello Rambaldi and F. Lillo.  The role of volume in order book dynamics: a multivariate Hawkes process analysis.  Quantitative Finance, 17(7):999–1020, 2017.
  * [27] H. Mei and J. M. Eisner.  The neural Hawkes process: a neurally self-modulating multivariate point process.  Advances in Neural Information Processing Systems, 30, 2017.
  * [28] G. O. Mohler, M. B. Short, P. J. Brantingham, F. P. Schoenberg, and G. E. Tita.  Self-exciting point process modeling of crime.  Journal of the American Statistical Association, 106(493):100–108, 2011.
  * [29] Y. Ogata.  The asymptotic behaviour of maximum likelihood estimators for stationary point processes.  Annals of the Institute of Statistical Mathematics, 30(1):243–261, 1978.
  * [30] Y. Ogata.  On lewis’ simulation method for point processes.  IEEE transactions on information theory, 27(1):23–31, 1981.
  * [31] Y. Ogata.  Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association, 83(401):9–27, 1988.
  * [32] Y. Ogata.  Space-time point-process models for earthquake occurrences.  Annals of the Institute of Statistical Mathematics, 50:379–402, 1998.
  * [33] Y. Ogata and H. Akaike.  On linear intensity models for mixed doubly stochastic Poisson and self- exciting point processes.  Journal of the Royal Statistical Society. Series B (Methodological), 44(1):102–107, 1982.
  * [34] Y. Ogata and J. Zhuang.  Space–time ETAS models and an improved extension.  Tectonophysics, 413(1):13–23, 2006.
  * [35] T. Omi, K. Aihara, et al.  Fully neural network based model for general temporal point processes.  Advances in Neural Information Processing Systems, 32, 2019.
  * [36] T. Ozaki.  Maximum likelihood estimation of Hawkes’ self-exciting point processes.  Annals of the Institute of Statistical Mathematics, 31:145–155, 1979.
  * [37] A. Reinhart.  A review of self-exciting spatio-temporal point processes and their applications.  Statistical Science, 33(3):299–318, 2018.
  * [38] M.-A. Rizoiu, Y. Lee, S. Mishra, and L. Xie.  Hawkes processes for events in social media.  In Frontiers of Multimedia Research, page 191–218. Association for Computing Machinery and Morgan & Claypool, 2017.
  * [39] M.-A. Rizoiu, S. Mishra, Q. Kong, M. Carman, and L. Xie.  SIR-Hawkes: Linking epidemic models and Hawkes processes to model diffusions in finite populations.  In Proceedings of the 2018 world wide web conference, pages 419–428, 2018.
  * [40] M.-A. Rizoiu, L. Xie, S. Sanner, M. Cebrian, H. Yu, and P. Van Hentenryck.  Expecting to be hip: Hawkes intensity processes for social media popularity.  In Proceedings of the 26th international conference on world wide web, pages 735–744, 2017.
  * [41] A. M. Schäfer and H. G. Zimmermann.  Recurrent neural networks are universal approximators.  In Artificial Neural Networks–ICANN 2006: 16th International Conference, Athens, Greece, September 10-14, 2006. Proceedings, Part I 16, pages 632–640. Springer, 2006.
  * [42] F. P. Schoenberg.  Multidimensional residual analysis of point process models for earthquake occurrences.  J. Amer. Stat. Assoc., 98(464):789–795, 2003.
  * [43] F. P. Schoenberg.  Testing separability in spatial-temporal marked point processes.  Biometrics, pages 471–481, 2004.
  * [44] F. P. Schoenberg.  Nonparametric estimation of variable productivity Hawkes processes.  Environmetrics, 33(6):e2747, 2022.
  * [45] F. P. Schoenberg, M. Hoffmann, and R. J. Harrigan.  A recursive point process model for infectious diseases.  Annals of the Institute of Statistical Mathematics, 71:1271–1287, 2019.
  * [46] O. Shchur, A. C. Türkmen, T. Januschowski, and S. Günnemann.  Neural temporal point processes: A review.  arXiv preprint arXiv:2104.03528, 2021.
  * [47] A. Veen and F. P. Schoenberg.  Estimation of space–time branching process models in seismology using an EM–type algorithm.  Journal of the American Statistical Association, 103(482):614–624, 2008.
  * [48] D. Vere-Jones.  Probabilities and information gain for earthquake forecasting.  Comput. Seismol., 30:248–263, 1998.
  * [49] D. Vere-Jones and T. Ozaki.  Some examples of statistical estimation applied to earthquake data.  Annals of the Institute of Statistical Mathematics, 34(1):189–207, 1982.
  * [50] T. Wang, M. Bebbington, and D. Harte.  Markov-modulated Hawkes process with stepwise decay.  Annals of the Institute of Statistical Mathematics, 64(3):521–544, 2012.
  * [51] L. Wasserman.  All of Nonparametric Statistics.  Springer New York, NY, 2006.
  * [52] J. Zhuang.  Next-day earthquake forecasts for the Japan region generated by the ETAS model.  Earth, Planets and Space, 63:207–216, 2011.
  * [53] J. Zhuang, T. Wang, and K. Kiyosugi.  Detection and replenishment of missing data in marked point processes.  Statistica Sinica, 30(4):2105–2130, 2020.
  * [54] J. Zhuang, M. J. Werner, and D. S. Harte.  Stability of earthquake clustering models: Criticality and branching ratios.  Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 88(6):062109, 2013.
  * [55] S. Zuo, H. Jiang, Z. Li, T. Zhao, and H. Zha.  Transformer Hawkes process.  In International Conference on Machine Learning, pages 11692–11702. PMLR, 2020.