Then, for any $M_T \to \infty$, $C > 0$ such that $\max_{l,k} \|h_{lk}^0\|_\infty \le C$, and $\hat{Q}$ defined in (9), we have that

$$ \mathbb{E}_0 \left[ \Pi \left( \|f^C - f_0\|_1 > M_T \sqrt{\kappa_T \epsilon_T} \right) \right] \xrightarrow{T \to \infty} 0 $$

$$ \mathbb{E}_0 \left[ \hat{Q} \left( \|f^C - f_0\|_1 > M_T \sqrt{\kappa_T \epsilon_T} \right) \right] \xrightarrow{T \to \infty} 0 $$

$$ \mathbb{P}_0 \left( \sum_{l,k=1}^{K} \left| \mathbb{E}_{\hat{Q}} [ \|h_{lk}^C\|_1 ] - \|h_{lk}^0\|_1 \right| > M_T \sqrt{\kappa_T \epsilon_T} \right) \xrightarrow{T \to \infty} 0. $$

**Proof** For the first statement on the posterior distribution, the argument of this proof is the same as the proof of Theorem 5.5 in Sulem et al. (2024) for the concentration in stochastic distance $\tilde{d}_{1T}$ and the following inequality obtained in the proof of Theorem 7 for the sigmoid link function:

$$ \mathbb{E}_0 \left[ \mathbb{1}_{\Omega_T} \Pi[f \in \mathcal{F} : \tilde{d}_{1T}(f, f_0) > M_T' \epsilon_T | N] \right] \le 4(2K + 1)e^{-x_1 M_T'^2 T \epsilon_T^2 / 2}, $$

where $x_1 > 0$, $M_T' > 0$, $\epsilon_T = \sqrt{\kappa_T \epsilon_T}$ and $\Omega_T$ is defined in the proof of Theorem 7. We now state and prove a result similar to Lemma 18, which holds if the inverse of the link function is Lipschitz on a large enough interval, for the specific case of the sigmoid link function.

**Lemma 23** Let $f \in \mathcal{F}_T$ such that $\|v - v_0\|_1 \le C$. Then

$$ \mathbb{E}_f [Z_{1l}] \ge \tilde{\zeta} \|f^C - f_0\|_1, $$

with $\tilde{\zeta} > 0$ a constant that may depend on $f_0$, C and $(\alpha, \eta, (\theta_k)_k)$.

The proof of this lemma follows the same argument as that of Lemma 4 in Sulem et al. (2024) (Section S9.20). With $x > 0$, for each $k \in [K]$, we define the event $\Omega_k$ as

$$ \Omega_k = \left\{ \max_{k' \neq k} N^{k'}[\tau_1, \tau_2) = 0, N^k[\tau_1, \tau_1 + x] = 0, N^k[\tau_1 + x, \tau_1 + x + A] = 1, N^k[\tau_1 + x + A, \tau_2) = 0 \right\}. $$

On $\Omega_k$, we have

$$ \mathbb{E}_f [Z_{1l}] \ge \sum_k \mathbb{E}_f \left[ \mathbb{1}_{\Omega_k} \int_{\tau_1}^{A+U_1^{(1)}} |\lambda_l^l(f) - \lambda_l^l(f_0)| dt \right], $$

where $U_1^{(1)}$ is the time of the first event on $N^k$. Using $\mathcal{Q}$ the point process measure of a homogeneous Poisson process with unit intensity on $\mathbb{R}^+$ and equal to the null measure on $[-A, 0)$, we can rewrite the previous inequality as

$$ \mathbb{E}_f [Z_{1l}] \ge \sum_k \mathbb{E}_\mathcal{Q} \left[ \int_{\tau_1}^{U_1^{(1)}+A} \mathcal{L}_l(f) \mathbb{1}_{\Omega_k} |\lambda_l^l(f) - \lambda_l^l(f_0)| dt \right], $$

with $\mathcal{L}_l(f)$ the likelihood process given by

$$ \mathcal{L}_l(f) = \exp \left( Kt - \sum_k \int_{\tau_1}^t \lambda_u^k(f) du + \sum_k \int_{\tau_1}^t \log(\lambda_u^k(f)) dN_u^k \right), \quad t \ge 0. $$
On $\Omega_k$, for $t \in [\tau_1, U_1^{(1)} + A)$, we have

$$
\begin{aligned}
\mathcal{L}_t(f) &\ge e^{Kt} \lambda_{U_1^{(1)}}^k(f) \exp \left\{ -\sum_{k'} \int_{\tau_1}^t \phi_{k'}(\tilde{\lambda}_{u}^{k'}(f)) du \right\} \\
&\ge \phi_k(v_k) \exp \left\{ -(x+A) \sum_k \theta_k \right\} \\
&\ge \phi_k(v_k^0 - C) \exp \left\{ -(x+A) \sum_k \theta_k \right\} = C_1,
\end{aligned}
$$

using that $f$ is such that $\|v - v_0\|_1 \le C$. Hence we have

$$
\mathbb{E}_f [Z_{1|1}] \ge C_1 \sum_k \mathbb{E}_{\mathbb{Q}} \left[ 1_{\Omega_k} \int_{U_1^{(1)}}^{U_1^{(1)}+A} \left| \phi_l(\tilde{\lambda}_t^l(f)) - \phi_l(\tilde{\lambda}_t^l(f_0)) \right| dt \right].
$$

Moreover, since for any $t \in [U_1^{(1)}, U_1^{(1)} + A]$, $|h_{kl}(t - U_1^{(1)} - h_{kl}^0(t - U_1^{(1)}))| \ge |h_{kl}^C(t - U_1^{(1)} - h_{kl}^0(t - U_1^{(1)}橋)|$,
we have

$$
\begin{aligned}
|\phi_l(\tilde{\lambda}_t^l(f)) - \phi_l(\tilde{\lambda}_t^l(f_0))| &= |\phi_l(v_l + h_{kl}(t - U_1^{(1)})) - \phi_l(v_l^0 + h_{kl}^0(t - U_1^{(1)}))| \\
&\ge |\phi_l(v_l + h_{kl}^C(t - U_1^{(1)})) - \phi_l(v_l^0 + h_{kl}^0(t - U_1^{(1)}))| = |\phi_l(\tilde{\lambda}_t^l(f^C)) - \phi_l(\tilde{\lambda}_t^l(f_0))|,
\end{aligned}
$$

and $\tilde{\lambda}_t^l(f^C) \in [-3C + 3C]$. On $[-3C + 3C]$, $\phi_l^{-1}$ is $L_{l,C}$-Lipschitz with $L_{l,C} > 0$, which leads to

$$
\begin{aligned}
\mathbb{E}_f [Z_{1|1}] &\ge C_1 \sum_k \mathbb{E}_{\mathbb{Q}} \left[ 1_{\Omega_k} \int_{U_1^{(1)}}^{U_1^{(1)}+A} \left| \phi_l(\tilde{\lambda}_t^l(f^C)) - \phi_l(\tilde{\lambda}_t^l(f_0)) \right| dt \right] \\
&\ge \frac{C_1}{L_{l,C}} \sum_k \mathbb{E}_{\mathbb{Q}} \left[ 1_{\Omega_k} \int_{U_1^{(1)}}^{U_1^{(1)}+A} \left| v_l + h_{kl}^C(t - U_1^{(1)} - v_l^0 - h_{kl}^0(t - U_1^{(1)})) \right| dt \right].
\end{aligned}
$$

Then the rest of proof is identical to that of Lemma A.4 in Sulem et al. (2024), replacing $h$ by $h^C$
and the constant $\tilde{\zeta}$ in the statement of the lemma can depend on $f_0$, $C$ and $L_{l,C}$.

We now prove Theorem 22. Let $f \in \mathcal{F}_T$ such that $\|f^C - f_0\|_1 > M_T \epsilon_T$ with $M_T$ such that $M_T \ge 4\mathbb{E}_0[\Delta\tau_1] M_T'$. Note that $\|f - f_0\|_1 \ge \|f^C - f_0\|_1$. Using Lemma 23 we have

$$
\begin{aligned}
\mathbb{E}_f [1_{\Omega_T} 1_{B_{d_1}(M_T' \epsilon_T)} | \mathcal{G}_0] &\le \mathbb{P}_f \left[ \sum_{j=1}^{J_T-1} Z_{jl} \le T M_T' \epsilon_T | \mathcal{G}_0 \right] \\
&\le \sum_{J \in \mathcal{J}_T} \mathbb{P}_f \left[ \sum_{j=1}^{J-1} Z_{jl} - \mathbb{E}_f [Z_{jl}] \le T M_T' \epsilon_T - \frac{T}{2 \mathbb{E}_0 [\Delta \tau_1]} \tilde{\zeta} M_T \epsilon_T | \mathcal{G}_0 \right] \\
&\le \sum_{J \in \mathcal{J}_T} \mathbb{P}_f \left[ \sum_{j=1}^{J-1} Z_{jl} - \mathbb{E}_f [Z_{jl}] \le - \frac{T}{4 \mathbb{E}_0 [\Delta \tau_1]} \tilde{\zeta} M_T \epsilon_T | \mathcal{G}_0 \right],
\end{aligned}
$$
since $M_T \ge 4\mathbb{E}_0[\Delta\tau_1] M'_T$. Besides, denoting $S_k = \{x : h_{kl}(x) \in [-2C, 2C]^c\}$ and $S_+ = \cup S_k$, we have

$$
\begin{aligned}
Z_{1l} &= \int_{\tau_1}^{\xi_1} |\phi(\tilde{\lambda}_l^t(f)) - \phi(\tilde{\lambda}_l^t(f_0))| dt = (U_1 - \tau_1)|\phi(v_l) - \phi(v_l^0| + \int_{U_1}^{\xi_1} |\phi(\tilde{\lambda}_l^t(f)) - \phi(\tilde{\lambda}_l^t(f_0))| dt \\
&\le (U_1 - \tau_1)|\phi(v_l) - \phi(v_l^0| + |S_+|\theta_l + \int_{U_1}^{U_1+A} |\phi(v_l + \sum_k 1_{U_1 \in N^k} h_{kl}^C(t - U_1)) - \phi(v_l^0| + \sum_k 1_{U_1 \in N^k} h_{kl}^0(t - U_1))| dt \\
&\le (U_1 - \tau_1)|\phi(v_l) - \phi(v_l^0| + |S_+|\theta_l + \int_0^A |\phi(v_l + \sum_k 1_{U_1 \in N^k} h_{kl}^C(t)) - \phi(v_l^0| + \sum_k 1_{U_1 \in N^k} h_{kl}^0(t))| dt.
\end{aligned}
$$

Moreover, since $\phi_l$ is $L_l$-Lipschitz with $L_l = \theta_l \alpha$, on the one hand we have

$$
\begin{aligned}
& \int_0^A |\phi(v_l + \sum_k 1_{U_1 \in N^k} h_{kl}^C(t) - \phi(v_l^0 + \sum_k 1_{U_1 \in N^k} h_{kl}^0(t))| dt \\
& \le L_l \int_0^A |v_l + \sum_k 1_{U_1 \in N^k} h_{kl}^C(t) - v_l^0 + \sum_k 1_{U_1 \in N^k} h_{kl}^0(t)| dt \\
& \le L_l (|v_l - v_l^0| + \sum_k \|h_{kl}^C - h_{kl}^0\|_1) \le L_l \|f^C - f_0\|_1.
\end{aligned}
$$

On the other hand we have $\|h_{kl}^C - h_{kl}^0\|_1 \ge |S_k|C$ which implies

$$
\|f^C - f_0\|_1 \ge \sum_k \|h_{kl}^C - h_{kl}^0\|_1 \ge \sum_k |S_k|C \ge C|S^+|.
$$

Hence we obtain

$$
Z_{1l} \le (U_1 - \tau_1) |\phi(v_l) - \phi(v_l^0)| + (\theta_l C^{-1} + L_l) \|f^C - f_0\|_1 \le [(U_1 - \tau_1)L_l + (\theta_l C^{-1} + L_l)] \|f^C - f_0\|_1.
$$

Therefore, since $\mathbb{E}_f[(\xi_1 - \tau_1)^n] = \frac{n!}{\|\phi(v)\|_1^n} \le \frac{n!}{\|\phi(v_0-C)\|_1^n}$, we obtain

$$
\begin{aligned}
\mathbb{E}_f [Z_{1l}^n] &\le 2^{n-1} [L_l^n \mathbb{E}_f [(\xi_1 - \tau_1)^n] + (\theta_l C^{-1} + L_l)^n] \|f^C - f_0\|_1^n \\
&\le 2^{n-1} [L_l^n \frac{n!}{\|\phi(v_0-C)\|_1^n} + (\theta_l C^{-1} + L_l)^n] \|f^C - f_0\|_1^n \\
&\le 2^{n-1} n! \max \left( \frac{L_1}{\|\phi(v_0-C)\|_1^n}, (\theta_l C^{-1} + L_l) \right)^n \|f^C - f_0\|_1^n \\
&\le \frac{1}{2} n! b^{n-2} v^2,
\end{aligned}
$$

with $b := 2 \max(\frac{L_1}{\|\phi(v_0-C)\|_1^n}, (\theta_l C^{-1} + L_l)) \|f^C - f_0\|_1$ and $v := 2 \max(\frac{L_1}{\|\phi(v_0-C)\|_1^n}, (\theta_l C^{-1} + L_l)) \|f^C - f_0\|_1$. Applying Bernstein's inequality as in Sulem et al. (2024), we obtain for $J \ge \frac{T}{2\mathbb{E}_0[\Delta\tau_1]}$,

$$
\mathbb{P}_f \left[ \sum_{j=1}^{J-1} Z_{jl} - \mathbb{E}_f [Z_{jl}] \le - \frac{T}{4\mathbb{E}_0[\Delta\tau_1]} \tilde{\zeta} M_T \varepsilon_T |\mathcal{G}_0| \right] \le \exp \left\{ - \frac{\tilde{\zeta}^2 T}{16 \bar{m}'} \right\}.
$$
since with $\bar{m} := \max\left(\frac{L_1}{\|\phi(v_0-C)\|_1^n}, (\theta_l C^{-1} + L_l)\right)$ we have

$$v^2 + b \frac{T}{4\mathbb{E}_0[\Delta\tau_1]} \bar{\zeta} M_T \varepsilon_T = 4\bar{m}^2 \|f^C - f_0\|_1^2 + \bar{m} \frac{T}{4\mathbb{E}_0[\Delta\tau_1]} \bar{\zeta} M_T \varepsilon_T \|f^C - f_0\|_1 \ge \bar{m}' M_T^2 \varepsilon_T^2.$$

and with $\bar{m}' = 4(\bar{m}^2 \vee \bar{m} \frac{T}{4\mathbb{E}_0[\Delta\tau_1]}) \bar{\zeta}$. Following the arguments of Sulem et al. (2024), this proves the first statement of Theorem 22 and the two other statements can be proved in the same way as Theorem 7 (see Appendix D.2), replacing $f$ by $f^C$ since in this case $\|h_{lk^C}\|_1 \le 2CA$ which implies that $\|f^C - f_0\|_1$ is bounded. ■

## Appendix F. Gibbs Sampler in the Sigmoid Hawkes Model

In this section, we describe a non-adaptive Gibbs sampler that computes the posterior distribution in the sigmoid Hawkes model, using the data augmentation scheme of Section 4 (see also Remark 5).

<table><caption>**Algorithm 4:** Gibbs sampler in the sigmoid Hawkes model with data augmentation</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>$N = (N^1, \dots, N^K)$, $n_{iter}$, $\mu$, $\Sigma$.</td></tr><tr><td>Output:</td><td>Samples $S = (f_i)_{i \in [n_{iter}]}$ from the posterior distribution $\Pi_A(f|N)$.</td></tr><tr><td>1</td><td>Precompute $(H_k(T_i^k))_{i \in [K]}$.</td></tr><tr><td>2</td><td>Initialize $f \sim N(f, \mu, \Sigma)$ and $S = []$.</td></tr><tr><td>3</td><td>for $t \leftarrow 1$ to $n_{iter}$ do</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;for $k \leftarrow 1$ to $K$ do</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for $i \leftarrow 1$ to $N_k$ do</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $\omega_i^k \sim p_{PG}(\omega_i^k; 1, \tilde{\lambda}_{T_i^k}^k(f))$</td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $(\tilde{T}_j^k)_{j=1,R_k}$ a Poisson temporal point process on $[0, T]$ with intensity $\theta_k \sigma(-\tilde{\lambda}_i^k(f))$</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;for $j \leftarrow 1$ to $R_k$ do</td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $\bar{\omega}_j^k \sim p_{PG}(\omega; 1, \tilde{\lambda}_{\bar{T}_j^k}^k(f))$</td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update $\tilde{\Sigma}_k = [\beta^2 H_k D_k (H_k)^T + \Sigma^{-1}]^{-1}$</td></tr><tr><td>11</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update $\tilde{\mu}_k = \tilde{\Sigma}_k (H_k [\beta v_k + \beta^2 \eta u_k] + \Sigma^{-1}\mu)$</td></tr><tr><td>12</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $f_k \sim N(f_k; \tilde{\mu}_k, \tilde{\Sigma}_k)$</td></tr><tr><td>13</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Add $f = (f_k)_k$ to $S$.</td></tr></tbody></table>

## Appendix G. Additional Results from our Numerical Experiments

In this section, we report results from our simulation study in Section 6 that were not added to the main text for conciseness purposes. Each of the following sub-sections corresponds to one of the simulation set-up.
<figure><img src="image_27.png" /><figcaption>Figure 25: Posterior distribution on $f = (v_1, h_{11})$ obtained with the MH sampler in the sigmoid model, in the three scenarios of Simulation 1 ($K = 1$). The three columns correspond to the *Excitation only* (left), *Mixed effect* (center), and *Inhibition only* (right) scenarios. The first row contains the marginal distribution on the background rate $v_1$, and the second row represents the posterior mean (solid orange line) and 95% credible sets (orange areas) on the (self) interaction function $h_{11}$. The true parameter $f_0$ is plotted in dotted green line.</figcaption></figure>

## G.1 Simulation 1

This section contains our results for the MH sampler, in the univariate settings of Simulation 1 with sigmoid and softplus link functions (see Figures 25 and 26).

## G.2 Simulation 3

This section contains additional results from Simulation 3: the estimated intensity function in the univariate and bivariate, well-specified settings (Figures 27 and 28), the estimated parameter in the mis-specified settings (Figure 29), and the estimated interaction functions in the bivariate setting and Inhibition scenario (Figure 30).
<figure><img src="image_28.png" /><figcaption>Figure 26: Posterior distribution on $f = (v_1, h_{11})$ obtained with the MH sampler in the softplus model, in the three scenarios of Simulation 1 ($K = 1$). The three columns correspond to the *Excitation only* (left), *Mixed effect* (center), and *Inhibition only* (right) scenarios. The first row contains the marginal distribution on the background rate $v_1$, and the second row represents the posterior mean (solid orange line) and 95% credible sets (orange areas) on the (self) interaction function $h_{11}$. The true parameter $f_0$ is plotted in dotted green line.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 27: Intensity function on a subwindow of the observation window estimated via the variational posterior mean and via the posterior mean computed with the MH sampler, in the well-specified setting of Simulation 3 on [0, 10], using the fully-adaptive mean-field variational (FA-MF-VI) algorithm (Algorithm 2). The true intensity $\lambda_t^1(f_0)$ is plotted in dotted green line.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 28: Estimated intensity function based on the (variational) posterior mean, in the well-specified and bivariate setting of Simulation 3 on [0, 10], using the fully-adaptive mean-field variational (FA-MF-VI) algorithm (Algorithm 2). The true intensity $\lambda_t(f_0)$ is plotted in dotted green line.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 29: Model-selection variational posterior distributions on $f = (v_1, h_{11})$ in the univariate sigmoid model and mis-specified setting of Simulation 3, evaluated by the fully-adaptive mean-field variational (FA-MF-VI) algorithm (Algorithm 2). The two columns correspond to a (mostly) *Excitation* (left) and a (mostly) *Inhibition* (right) settings. The first row contains the marginal distribution on the background rate $v_1$, and the second row represents the variational posterior mean (solid line) and 95% credible sets (colored areas) on the (self) interaction function $h_{11}$. The true parameter $f_0$ is plotted in dotted green line.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 30: Posterior and model-selection variational posterior distributions on $f = (v, h)$ in the bivariate sigmoid model, well-specified setting, and Inhibition setting of Simulation 3, evaluated by the non-adaptive MH sampler and the fully-adaptive mean-field variational (FA-MF-VI) algorithm (Algorithm 2). The first row contains the marginal distribution on the background rates ($v_1, v_2$), and the second and third rows represent the (variational) posterior mean (solid line) and 95% credible sets (colored areas) on the four interaction function $h_{11}, h_{12}, h_{21}, h_{22}$. The true parameter $f_0$ is plotted in dotted green line.</figcaption></figure>
G.3 Simulation 4

This section contains our results for the Inhibition setting of Simulation 4, i.e., the estimated graphs
in (Figures 31 and 32), the heatmaps of the risk on the interaction functions in Figure 33, the
estimated L₁-norms after the first step of Algorithm 3 in Figure 34, and the variational posterior
distribution on the subset of the parameter in Figure 35. We also report a comparison of the risks
obtained after the first and second steps of Algorithm 3 in Figure 36.
<figure><img src="image_33.png" /><figcaption>Figure 31: Estimated graph parameter $\hat{\delta}$ (black=0, white=1) for $K = 2, 4, 8, 16, 32, 64$ in the Excitation scenario of Simulation 4.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure 32: Estimated graph parameter $\hat{\delta}$ (black=0, white=1) for $K = 2, 4, 8, 16, 32, 64$ in the Inhibition scenario of Simulation 4.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 33: Heatmaps of the $L_1$-norms of the true $h_0$, i.e., the entries of the matrix $S_0 = (S_{lk}^0)_{l,k} = (\|h_{lk}^0\|_1)_{l,k}$ (left column) and the $L_1$-error of the model selection variational posterior mean (right column) obtained with Algorithm 3, in the Inhibition scenario of Simulation 4. The rows correspond to $K = 2, 4, 8, 16, 32, 64$.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Figure 34: Estimated $L_1$-norms after the first step of Algorithm 3 (in blue), and ground-truth norms (in orange), plotted in increasing order, in the Inhibition scenario of Simulation 4, for the models with $K \in \{2, 4, 8, 16, 32, 64\}$.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Figure 35: Model-selection variational posterior distributions on $v_1$ (left column) and interaction functions $h_{11}$ and $h_{21}$ (second and third columns) in the Inhibition scenario and multi-variate sigmoid models of Simulation 4, computed with our two-step mean-field variational (MF-VI) algorithm (Algorithm 3). The different rows correspond to different multivariate settings $K = 2, 4, 8, 16, 32, 64$.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Figure 36: Risks of the variational posteriors obtained after the first step (in blue) and the second step (in orange) of Algorithm 3 in the *Excitation* scenario (left) and *Inhibition* scenario (right) of Simulation 4, for the models with $K = 2, 4, 8, 16, 32$. For each setting, we repeat the experiment 10 times and we plot the mean risk and the intervals at $\pm 2$ standard deviations.</figcaption></figure>
G.4 Simulation 5

In this section, we report some characteristics of the simulated data in Simulation 5, in particular
the number of points and excursions in each setting (see Table 9). Moreover, we report the plots of
the posterior distribution in a subset of the parameter in Figure 37.
<table><caption>Table 9: Number of points and *global* and average *local* excursions in the multidimensional data sets of Simulation 5 ($K = 10$).</caption><thead><tr><th>Scenario</th><th>T</th><th># events</th><th># excursions</th><th># local excursions</th></tr></thead><tbody><tr><td rowspan="4">Excitation</td><td>50</td><td>2621</td><td>36</td><td>114</td></tr><tr><td>200</td><td>10,729</td><td>155</td><td>473</td></tr><tr><td>400</td><td>21,727</td><td>303</td><td>957</td></tr><tr><td>800</td><td>42,904</td><td>596</td><td>1921</td></tr><tr><td rowspan="4">Inhibition</td><td>50</td><td>1747</td><td>49</td><td>134</td></tr><tr><td>200</td><td>7019</td><td>222</td><td>529</td></tr><tr><td>400</td><td>13,819</td><td>466</td><td>1053</td></tr><tr><td>800</td><td>27,723</td><td>926</td><td>2118</td></tr></tbody></table>

<figure><img src="image_39.png" /><figcaption>Figure 37: Model-selection variational posterior on two interaction functions $h_{66}$ and $h_{76}$, for different observation lengths $T \in \{50, 200, 400, 800\}$, in the *Excitation* and *Inhibition* scenarios in Simulation 5 with $K = 10$. We note that in this simulation, the true number of basis functions is 2 and is well recovered for all values of $T$. The estimation of these two interaction functions is poor for the smallest $T$, however, it improves when $T$ increases.</figcaption></figure>
G.5 Simulation 6

This section contains the estimated graphs (Figures 38 and 40), the variational posterior distribution
on a subset of the parameter (Figures 39 and 41), in the mis-specified settings of Simulation 6.
<figure><img src="image_40.png" /><figcaption>Figure 38: Estimated graph after thresholding the L₁-norms using the “gap” or “slope change” heuristic, in the different settings of mis-specified link functions of Simulation 6, and in the *Excitation* and *Inhibition* scenarios. We observe that the true graph (with non-null principal and first off-diagonal) is correctly estimated for the ReLU mis-specification setting, while some errors happen in the two other link settings, in particular in the *Inhibition* scenario.</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Figure 39: Estimated interaction functions $h_{66}$ and $h_{76}$ in the mis-specified settings of Simulation 6, where the data is generated from a Hawkes model with ReLU, softplus, or a mis-specified link function, and in the *Excitation* and *Inhibition* scenarios. We note that the estimation of the interaction functions is deteriorated in these mis-specified cases, however the sign of the functions are still recovered.</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>Figure 40: Estimated graph after thresholding the L₁-norms, when using Algorithm 3 with different support upper bounds $A' \in \{0.5, 0.1, 0.2, 0.4\}$, containing the true memory parameter $A = 0.1$, in the settings of Simulation 7. We note that the true graph (with non-null principal and first off-diagonal) is correctly estimated in all cases, in the *Excitation* scenario (first row) and in the *Inhibition* scenario (second row).</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>Figure 41: Estimated background rates $v_k$ for $k = 1, \dots, 5$ when using different values of the upper bound parameter $A \in \{0.05, 0.1, 0.2, 0.4\}$, in the two scenarios of Simulation 8. As expected, the background rates are better estimated in the well-specified setting $A = A_0 = 0.1$; nonetheless, when $A$ is not too far above $A_0$, the estimation does not deteriorate too much, in particular in the *Inhibition* scenarios.</figcaption></figure>

## References

Ryan Prescott Adams, Iain Murray, and David J. C. MacKay. Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities. In *Proceedings of the 26th Annual International Conference on Machine Learning*, ICML '09, page 9–16, 2009. doi: 10.1145/1553374.1553376.

Julyan Arbel, Ghislaine Gayraud, and Judith Rousseau. Bayesian optimal adaptive estimation using a sieve prior. *Scandinavian journal of statistics*, 40(3):549–570, 2013.

Emmanuel Bacry and Jean-François Muzy. Second order statistics characterization of Hawkes processes and non-parametric estimation, 2015.

Emmanuel Bacry, Martin Bompaire, Stéphane Gaiffas, and Jean-François Muzy. Sparse and low-rank multivariate Hawkes processes. *Journal of Machine Learning Research*, 21(50):1–32, 2020.

Christopher M. Bishop. *Pattern recognition and machine learning*. Information Science and Statistics. Springer, New York, 2006. doi: 10.1007/978-0-387-45528-0.

David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. *Journal of the American Statistical Association*, 112(518):859–877, Apr 2017. doi: 10.1080/01621459.2017.1285773.
Anna Bonnet, Miguel Martinez Herrera, and Maxime Sangnier. Maximum likelihood estimation for Hawkes processes with self-excitation or inhibition. *Statistics & Probability Letters*, 179:109214, 2021.

Pierre Brémaud and Laurent Massoulié. Stability of nonlinear Hawkes processes. *The Annals of Probability*, pages 1563–1588, 1996.

Biao Cai, Jingfei Zhang, and Yongtao Guan. Latent network structure learning from high-dimensional multivariate point processes. *Journal of the American Statistical Association*, 119 (545):95–108, 2024.

Chris U Carmona and Geoff K Nicholls. Scalable semi-modular inference with variational meta-posteriors. *arXiv preprint arXiv:2204.00296*, 2022.

Lisbeth Carstensen, Albin Sandelin, Ole Winther, and Niels R Hansen. Multivariate Hawkes process models of the occurrence of regulatory elements. *BMC bioinformatics*, 11(1):1–19, 2010.

I. Castillo and A. van der Vaart. Needles and straw in a haystack: Posterior concentration for possibly sparse sequences. *The Annals of Statistics*, 40:2069–2101, 2012.

I. Castillo, J. Schmidt-Hieber, and A. van der Vaart. Bayesian linear regression with sparse priors. *The Annals of Statistics*, 43:1986–2018, 2015.

Shizhe Chen, Ali Shojaie, Eric Shea-Brown, and Daniela Witten. The multivariate Hawkes process in high dimensions: Beyond mutual excitation. *arXiv:1707.04928v2*, 2017.

Manon Costa, Carl Graham, Laurence Marsalle, and Viet Chi Tran. Renewal in Hawkes processes with self-excitation and inhibition. *Advances in Applied Probability*, 52(3):879–915, 2020. doi: 10.1017/apr.2020.19.

Daryl J Daley and David Vere-Jones. *An introduction to the theory of point processes: volume II: general theory and structure*. Springer Science & Business Media, 2007.

Isabella Deutsch and Gordon J. Ross. Bayesian estimation of multivariate hawkes processes with inhibition and sparsity, 2022.

Christian Donner and Manfred Opper. Efficient bayesian inference of sigmoidal Gaussian Cox processes. *Journal of Machine Learning Research*, 19(67):1–34, 2018.

Sophie Donnet, Vincent Rivoirard, and Judith Rousseau. Nonparametric bayesian estimation for multivariate hawkes processes. *Annals of Statistics*, 48(5):2698–2727, 2020.

Michael Eichler, Rainer Dahlhaus, and Johannes Dueck. Graphical modeling for multivariate Hawkes processes with nonparametric link functions. *Journal of Time Series Analysis*, 38(2): 225–242, 2017.

Felipe Gerhard, Moritz Deger, and Wilson Truccolo. On the stability and dynamics of stochastic spiking neuron models: Nonlinear Hawkes process and point process glms. *PLOS Computational Biology*, 13:1–31, 02 2017. doi: 10.1371/journal.pcbi.1005390.
Gene H Golub and John H Welsch. Calculation of Gauss quadrature rules. *Mathematics of computation*, 23(106):221–230, 1969.

Niels Richard Hansen, Patricia Reynaud-Bouret, and Vincent Rivoirard. Lasso and probabilistic inequalities for multivariate point processes. *Bernoulli*, 21(1):83–143, 2015.

Alan G Hawkes. Point spectra of some mutually exciting point processes. *Journal of the Royal Statistical Society: Series B (Methodological)*, 33(3):438–443, 1971.

Alan G. Hawkes. Hawkes processes and their applications to finance: a review. *Quantitative Finance*, 18(2):193–198, 2018. doi: 10.1080/14697688.2017.1403131.

M. Hoffmann, J. Rousseau, and J. Schmidt-Hieber. On adaptive posterior concentration rates. *The Annals of Statistics*, 43:2259–2295, 2015.

J. F. C. Kingman. *Poisson processes*, volume 3 of *Oxford Studies in Probability*. The Clarendon Press Oxford University Press, New York, 1993.

Remi Lemonnier and Nicolas Vayatis. Nonparametric Markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate Hawkes processes. In *Joint European Conference on Machine Learning and Knowledge Discovery in Databases*, pages 161–176. Springer, 2014.

Rafael Lima. Hawkes processes modeling, inference, and control: An overview. *SIAM Review*, 65(2):331–374, 2023. doi: 10.1137/21M1396927.

Xiaofei Lu and Frédéric Abergel. High-dimensional Hawkes processes for limit order books: modelling, empirical analysis and numerical calibration. *Quantitative Finance*, 18(2):249–264, 2018.

Noa Malem-Shinitski, Cesar Ojeda, and Manfred Opper. Nonlinear Hawkes process with Gaussian process self effects, 2021.

Hongyuan Mei and Jason Eisner. The neural Hawkes process: A neurally self-modulating multivariate point process, 2017.

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. doi: 10.1198/jasa.2011.ap09546.

Dennis Nieman, Botond Szabo, and Harry van Zanten. Contraction rates for sparse variational approximations in Gaussian process regression. *Journal of Machine Learning Research*, 23(205):1–26, 2022.

Yosihiko Ogata. Seismicity analysis through point-process modeling: A review. *Seismicity patterns, their statistical significance and physical meaning*, pages 471–507, 1999.

Ilsang Ohn and Lizhen Lin. Adaptive variational Bayes: Optimality, computation and applications. *The Annals of Statistics*, 52(1):335–363, 2024.
Jack Olinde and Martin B. Short. A self-limiting Hawkes process: Interpretation, estimation, and use in crime modeling. In *2020 IEEE International Conference on Big Data (Big Data)*, pages 3212–3219, 2020. doi: 10.1109/BigData50022.2020.9378017.

Peter Pfaffelhuber, Stefan Rotter, and Jakob Stiefel. Mean-field limits for non-linear Hawkes processes with excitation and inhibition. *Stochastic Processes and their Applications*, 153:57–78, 2022.

Nicholas G Polson, James G Scott, and Jesse Windle. Bayesian inference for logistic models using pólya–gamma latent variables. *Journal of the American statistical Association*, 108(504):1339–1349, 2013.

Kolyan Ray and Botond Szabó. Variational Bayes for high-dimensional linear regression with sparse priors. *Journal of the American Statistical Association*, pages 1–12, jan 2021. doi: 10.1080/01621459.2020.1847121.

Weining Shen and Subhashis Ghosal. Adaptive Bayesian procedures using random series priors. *Scandinavian Journal of Statistics*, 42(4):1194–1213, 2015. doi: https://doi.org/10.1111/sjos.12159.

Deborah Sulem, Vincent Rivoirard, and Judith Rousseau. Bayesian estimation of nonlinear Hawkes processes. *Bernoulli*, 30(2):1257–1286, 2024.

Michalis Titsias and Miguel Lázaro-Gredilla. Spike and slab variational inference for multi-task and multiple kernel learning. In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K.Q. Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 24, 2011.

A. W. van der Vaart and J. H. van Zanten. Adaptive Bayesian estimation using a Gaussian random field with inverse Gamma bandwidth. *The Annals of Statistics*, 37(5B), oct 2009.

John Worrall, Raiha Browning, Paul Wu, and Kerrie Mengersen. Fifty years later: new directions in Hawkes processes. SORT (*Statistics and Operations Research Transactions*), 46(1):3–38, 2022.

Fengshuo Zhang and Chao Gao. Convergence rates of variational posterior distributions. *The Annals of Statistics*, 48(4):2180–2207, 2020.

Rui Zhang, Christian Walder, and Marian-Andrei Rizoiu. Variational inference for sparse Gaussian process modulated Hawkes process. *Proceedings of the AAAI Conference on Artificial Intelligence*, 34(04):6803–6810, Apr 2020. doi: 10.1609/aaai.v34i04.6160.

Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. Efficient inference for nonparametric Hawkes processes using auxiliary latent variables. *Journal of Machine Learning Research*, 21(241):1–31, 2020.

Feng Zhou, Quyu Kong, Yixuan Zhang, Cheng Feng, and Jun Zhu. Nonlinear Hawkes processes in time-varying system, 2021.

Feng Zhou, Quyu Kong, Zhijie Deng, Jichao Kan, Yixuan Zhang, Cheng Feng, and Jun Zhu. Efficient inference for dynamic flexible interactions of neural populations. *Journal of Machine Learning Research*, 23(211):1–49, 2022.