Goodness-of-Fit Test for Self-Exciting Processes

Recently there have been many research efforts in developing generative models for self-exciting point processes, partly due to their broad applicability for real-world applications, notably self- and mutual- exciting point processes. However, rarely can we quantify how well the generative model captures the nature or ground-truth since it is usually unknown. The challenge typically lies in the fact that the generative models typically provide, at most, good approximations to the ground-truth (e.g., through the rich representative power of neural networks), but they cannot be precisely the ground-truth. We thus cannot use the classic goodness-of-fit test framework to evaluate their performance. In this paper, we provide goodness-of-fit tests for generative models by leveraging a new connection of this problem with the classical statistical theory of mismatched maximum-likelihood estimator (MLE). We present a non-parametric self-normalizing test statistic for the goodness-of-fit test based on Generalized Score (GS) statistics. We further establish asymptotic properties for MLE of the Quasi-model (Quasi-MLE), asymptotic $\chi^2$ null distribution and power function of GS statistic. Numerical experiments validate the asymptotic null distribution as well as the consistency of our proposed GS test.


Introduction
Self-and mutual-exciting point processes, as known as the Hawkes processes, introduced by the original papers by Hawkes Hawkes [1971a,b], Hawkes and Oakes [1974], become popular in machine learning due to their wide applicability in modeling triggering effect in discrete event data is ubiquitous from modern applications ranging from seismology Ogata [1988,1999], Zhuang [2011], infectious disease modeling Meyer et al. [2014], Schoenberg et al. [2019], crime events Mohler et al. [2011], wildfire occurrence Peng et al. [2005], civilian deaths in Iraq Lewis et al. [2012], terrorist activity forecasting Porter et al. [2012], social network analysis and so on.
Classical Hawkes processes are largely parametric, which focus on modeling the conditional intensity function of the point process (since the conditional intensity function completely specifies the distribution of the process). Hawkes process assumes that the intensity function consists of the sum of a deterministic background intensity (which can be time-varying) and a stochastic term, which captures the influence from the past events. It is common to assume that the influence from past events is additive, and the influence of an individual event is measured by the so-called triggering function.
One key problem in the Hawkes process is to specify the triggering kernel. Popular parametric triggering functions include exponential kernel, power kernel, and Matérn kernel Reinhart et al. [2018]. g < l a t e x i t s h a 1 _ b a s e 6 4 = " D V Y S O L A d d / + O 0 f f + K Y T 5 C 6 L L R z k = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q Y I 8 F L z 1 W s B / Q h r L Z b t q l m 0 3 Y n Q g l 9 E d 4 8 a C I V 3 + P N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K F Q d f 9 d g p b 2 z u 7 e 8 X 9 0 s H h 0 f F J + f S s Y + J U M 9 5 m s Y x 1 L 6 C G S 6 F 4 G w V K 3 k s 0 p 1 E g e T e Y 3 i / 8 7 h P X R s T q E W c J 9 y M 6 V i I U j K K V u o M m x W w 8 H 5 Y r b t V d g m w S L y c V y N E a l r 8 G o 5 i l E V f I J D W m 7 7 k J + h n V K J j k 8 9 I g N T y h b E r H v G + p o h E 3 f r Y 8 d 0 6 u r D I i Y a x t K S R L 9 f d E R i N j Z l F g O y O K E 7 P u L c T / v H 6 K Y d 3 P h E p S 5 I q t F o W p J B i T x e 9 k J D R n K G e W U K a F v Z W w C d W U o U 2 o Z E P w 1 l / e J J 1 a 1 b u p 1 h 5 u K 4 1 6 H k c R L u A S r s G D O 2 h A E 1 r Q B g Z T e I Z X e H M S 5 8 V 5 d z 5 W r Q U n n z m H P 3 A + f w B j k 4 + S < / l a t e x i t >ĝ < l a t e x i t s h a 1 _ b a s e 6 4 = " D V Y S O L A d d / + O 0 f f + K Y T 5 C 6 L L R z k = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q Y I 8 F L z 1 W s B / Q h r L Z b t q l m 0 3 Y n Q g l 9 E d 4 8 a C I V 3 + P N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K F Q d f 9 d g p b 2 z u 7 e 8 X 9 0 s H h 0 f F J + f S s Y + J U M 9 5 m s Y x 1 L 6 C G S 6 F 4 G w V K 3 k s 0 p 1 E g e T e Y 3 i / 8 7 h P X R s T q E W c J 9 y M 6 V i I U j K K V u o M m x W w 8 H 5 Y r b t V d g m w S L y c V y N E a l r 8 G o 5 i l E V f I J D W m 7 7 k J + h n V K J j k 8 9 I g N T y h b E r H v G + p o h E 3 f r Y 8 d 0 6 u r D I i Y a x t K S R L 9 f d E R i N j Z l F g O y O K E 7 P u L c T / v H 6 K Y d 3 P h E p S 5 I q t F o W p J B i T x e 9 k J D R n K G e W U K a F v Z W w C d W U o U 2 o Z E P w 1 l / e J J 1 a 1 b u p 1 h 5 u K 4 1 6 H k c R L u A S r s G D O 2 h A E 1 r Q B g Z T e I Z X e H M S 5 8 V 5 d z 5 W r Q U n n z m H P 3 A + f w B j k 4 + S < / l a t e x i t > g ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 n h V t e 8 E 5 Q j H 0 7 E n q 3 M h u 3 O b d + c = " > A A A B 6 n i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R Z B P J T d K t h j w Y v H i v Y D 2 r V k 0 + w 2 N J s s S V Y o S 3 + C F w + K e P U X e f P f m L Z 7 0 N Y H A 4 / 3 Z p i Z F y S c a e O 6 3 0 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t b V M F a E t I r l U 3 Q B r y p m g L c M M p 9 1 E U R w H n H a C 8 c 3 M 7 z x R p Z k U D 2 a S U D / G k W A h I 9 h Y 6 T 5 6 v B i U K 2 7 V n Q O t E i 8 n F c j R H J S / + k N J 0 p g K Q z j W u u e 5 i f E z r A w j n E 5 L / V T T B J M x j m j P U o F j q v 1 s f u o U n V l l i E K p b A m D 5 u r v i Q z H W k / i w H b G 2 I z 0 s j c T / / N 6 q Q n r f s Z E k h o q y G J R m H J k J J r 9 j Y Z M U W L 4 x B J M F L O 3 I j L C C h N j 0 y n Z E L z l l 1 d J u 1 b 1 L q u 1 u 6 t K o 5 7 H U Y Q T O I V z 8 O A a G n A L T W g B g Q i e 4 R X e H O 6 8 O O / O x 6 K 1 4 O Q z x / A H z u c P 4 + a N g Q = = < / l a t e x i t > g 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " C e 9 5 s D 4 R / U g L i G s c 4 z r r B N w 3 F t Q = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q Y I 8 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p Y T R w B + W K W 3 U X I O v E y 0 k F c j Q H 5 a / + M G Z p h N I w Q b X u e W 5 i / I w q w 5 n A W a m f a k w o m 9 A R 9 i y V N E L t Z 4 t T Z + T C K k M S x s q W N G S h / p 7 I a K T 1 N A p s Z 0 T N W K 9 6 c / E / r 5 e a s O 5 n X C a p Q c m W i 8 J U E B O T + d 9 k y B U y I 6 a W U K a 4 v Z W w M V W U G Z t O y Y b g r b 6 8 T t q 1 q n d V r d 1 f V x r 1 P I 4 i n M E 5 X I I H N 9 C A O 2 h C C x i M 4 B l e 4 c 0 R z o v z 7 n w s W w t O P n M K f + B 8 / g D u g 4 2 I < / l a t e x i t > g ⇤ = g 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 O D 9 V C Y l + V r r o U 9 S / 9 G 8 z k v D w 7 w = " > A A A B 7 n i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I s g H s p u F e x F K H j x W M F + Q L u W b J r d h m a z S 5 I V y t I f 4 c W D I l 7 9 P d 7 8 N 6 b t H r T 1 w c D j v R l m 5 v m J 4 N o 4 z j c q r K 1 v b G 4 V t 0 s 7 u 3 v 7 B + X D o 7 a O U 0 V Z i 8 Y i V l 2 f a C a 4 Z C 3 D j W D d R D E S + Y J 1 / P H t z O 8 8 M a V 5 L B / M J G F e R E L J A 0 6 J s V I n f L y 4 C Q f O o F x x q s 4 c e J W 4 O a l A j u a g / N U f x j S N m D R U E K 1 7 r p M Y L y P K c C r Y t N R P N U s I H Z O Q 9 S y V J G L a y + b n T v G Z V Y Y 4 i J U t a f B c / T 2 R k U j r S e T b z o i Y k V 7 2 Z u J / X i 8 1 Q d 3 L u E x S w y R d L A p S g U 2 M Z 7 / j I V e M G j G x h F D F 7 a 2 Y j o g i 1 N i E S j Y E d / n l V d K u V d 3 L a u 3 + q t K o 5 3 E U 4 Q R O 4 R x c u I Y G 3 E E T W k B h D M / w C m 8 o Q S / o H X 0 s W g s o n z m G P 0 C f P 0 6 r j t w = < / l a t e x i t > Figure 1: The ground truth is g * ; the assumed family of candidate models is G. Goodness-of-fit addresses how close the fitted model g is to the unknown true one g * . In classic set up, one assumes there exists a g 0 ∈ G such that g * = g 0 . Here, we treat a more general case where g * may not be contained in G, but there exists a good approximate g 0 ∈ G to it. Under this setting, goodness-of-fit answers how close g is g 0 instead.
(2) There still lacks a principled method to quantify the goodness-of-fit for triggering effect in the Hawkes processes, which is the main effectof-interest in many Hawkes process models. Existing goodness-of-fit tests for Hawkes processes, such as Ogata [1988[ ], Schoenberg [2003 aim to test the whole conditional intensity function. Since the background rate is usually the dominating term in the conditional intensity, this test may not detect subtle triggering function differences.
In this paper, we present a non-parametric goodness-of-fit test statistic, called the Generalized Score (GS) statistic, which can be broadly applied to evaluating self-exciting process parametric models and generative models. The GS test is constructed by reformulating the goodness-of-fit test as a two-sample test: whether the real data and data generated from the generative model have the same distribution? Based on this, we derive the likelihood score statistic with estimated piece-wise constant kernels to fully adapt to data and with minimum model restrictions. We further establish asymptotic properties for MLE of the Quasi-model (Quasi-MLE), asymptotic χ 2 null distribution, and power function of GS statistic. The main ingredients of our analysis include (1) making a connection between GS test and the classic theory on MLE under model misspecification White [1982] and (2) generalizing the asymptotic properties of MLE of Hawkes process in Ogata et al. [1978] to model misspecification case. Our GS test provides a tool for model diagnosis and comparison of Hawkes process generative models. Using numerical simulation and real-data examples, we demonstrate the effectiveness of our proposed test.
Several features of our GS test include: (1) We develop the test for generative models considering their inherent "model misspecification nature"; (2) we focus on goodness-of-fit of the triggering effect in Hawkes process models; (3) due to its construction, the GS statistic enjoys simple asymptotic distribution specified by χ 2 distribution and analytical form of the power function, which may enable us to select a threshold for the test and calibrate the test.
Related Work. One-sample goodness-of-fit problem is closely related to the two-sample test problem. For independent and identically distributed (i.i.d.) observations, two-sample test is well studied (e.g. energy statistic Székely et al. [2004], Baringhaus and Franz [2004] and maximum mean discrepancy (MMD) Gretton et al. [2012]) and so is the goodness-of-fit based on it. Chwialkowski et al. [2016] developed Stein operator based MMD (which they call squared Stein discrepancy) and changed the two-sample test statistic to a one-sample goodness-of-fit test metric. Bounliphone et al. [2015] reformulated the one-sample goodness-of-fit problem into a two-sample test problem and developed a model selection tool based on MMD. Extension of those methods to point process is missing until Yang et al. [2019] proposed a kernel goodness-of-fit test by defining a Stein discrepancy for generic point process; However, a common drawback of a kernel-based test is that the null distribution is hard to evaluate (since they depend on infinite series involving the eigenvalues of the kernel). In contrast, our GS statistic follows a simple χ 2 null distribution and is easy to calibrate. Our proposed method allows the distribution under the null to be flexible and estimated from data by comparing the data to the generative model via the test statistic. Other model diagnostics include likelihood of fitted model and the observed data Schorlemmer et al. [2007] and Information Criterion (IC) Chen et al. [2018]. The likelihood is the most commonly used, but overfitting makes it less convincing and even questionable. Chen et al. [2018] assumed correct model specification, which typically does not hold in the real study, and the consistency result of IC is restricted to exponential triggering function case. For more on the kernel-based two-sample test as well as model diagnosis and selection method of the point process, one can refer to Harchaoui et al. [2013] and Bray and Schoenberg [2013].

Mathematical Background
Consider a counting process {N (t) : t ≥ 0}, with associated history H 0,t = {t i : 0 < t i < t} (t ≥ 0) indicating the occurrence time of a sequence of discrete events. We use H t for H 0,t for simplification. A point process is characterized by its conditional intensity function, which defined as follows: Hawkes process is a self-exciting point process with conditional intensity takes the following form: where µ is called the background intensity and φ : (0, ∞) → [0, ∞) is called the triggering function.
To simplify the notation, we can define φ on R but it takes value zero on (−∞, 0].
We assume the separability of triggering function into components for magnitude and time: where temporal triggering function g is a probability density function (p.d.f.) and α represents the magnitude of triggering effect, i.e. how many subsequent events one event can trigger on average. Given the past trajectory H T with N events, the log-likelihood over time interval [0, T ] can be expressed as: One can refer to Laub et al. [2015] and Reinhart et al. [2018] for a more comprehensive introduction of Hawks process as well as a detailed deviation of its (log-)likelihood function.

A Non-Parametric Goodness-of-Fit Test
In this section, we propose a diagnosis procedure for goodness-of-fit for fitted Hawkes process generative models based on Generalized Score test under model mis-specification. We first reformulate the one-sample goodness-of-fit problem into a two-sample testing problem.

Problem set-up
Suppose we have two data sequences D z = (t Nz ), (z = 1, 2), which represent the arrival time of a sequence of events. Here, D 1 is from real world and D 2 is generated from the fitted generative model. Assume D 1 ∼ λ * and D 2 ∼ λ, where λ * is the true but unknown conditional intensity and λ is the fitted one. Further assume both conditional intensities take form in (1). We aim to test H 0 : φ * = φ, versus H 1 : φ * = φ. Note that φ * in the above formulation is unknown. In the following of this section, we'll propose a two-sample test statistic that helps us to quantify the goodness-of-fit of φ. What's more, if we fit another generative model, we can compare those two models concerning their triggering components.
We calculate this test statistic in the following three steps: Mix the two data sequences up to get an aggregated sequence; Estimate θ 0 , maximizer of the Quasi-likelihood, from a Quasi-parameter space Θ for the aggregated sequence; Compute a test statistic GS T based on the estimation in the last step.
Remark 1 (Singleton null). In our setting, the unknown parameter of the triggering function is infinite-dimensional, so the null hypothesis is an uncountable set. To make the problem tractable, we simplify by using singleton null hypothesis, which corresponds to the two triggering function are identical (this is a typical way to set-up a two-sample test. Moreover, we will represent the unknown triggering function using some basis function (in our case, we use indicator function on mutually disjoint intervals, see (2)) such that we reduce this into a finite-dimensional problem.
Remark 2 (Model mismatch). We use the term "Quasi" here, since commonly speaking there will be a mismatch between a machine learning algorithm class we specify and the unknown true intensity, i.e., this class is misspecified as illustrated in Figure 1. We add a prefix "Quasi-" for everything under this class, e.g., Quasi-conditional intensity and Quasi-likelihood function. Since conditional intensity characterizes a point process and we assume the triggering function takes the form φ * = αg * , we only need to specify the approximate class for g * . Here we use piecewise constant function class G, since a piecewise constant function can be arbitrarily close to an integrable function when reducing the size of the discretization bin. Most importantly, empirically, there exists g 0 ∈ G, which corresponds to our estimand θ 0 , which serves as a good approximation to g * .

Proposed goodness-of-fit test procedure
The idea behind this test comes from a critical observation that under H 0 , mixing two sequences will lead to a Hawkes process with scaled intensity function. Based on this observation, we can derive a Generalized Score (GS) test, which is known to be locally most powerful (Neyman-Pearson lemma).
Step 1: mix two data sequences and model the aggregated sequence. In this step, we derive the Quasi-log-likelihood function for the aggregated sequence. The proof is in Appendix A. Proposition 1 (Log-likelihood of mixing of two Hawkes processes). Suppose we have two Hawkes processes with conditional intensities Define their mixing to be N (t) = N 1 (t) + N 2 (t). Then it has background intensity µ = µ (1) + µ (2) . Denote T = max{t Nz }, z = 1, 2, we have that: (i) Under H 1 , let z = 2(or 1) when z = 1(or 2), the full model log- Remark. (i) Note that the triggering function takes value zero on (−∞, 0] and thus we did not consider the triggering effect of events to its own history. (ii) Under H 0 , we can model the aggregated data using a univariate Hawkes process with the same excitation function. (iii) For each event in process z (z = 1, 2), it does not only dependent its original own history, but also depends on the history of another Hawkes process z (z = 2, 1). An illustration of this is given in Figure 5 in Appendix A.
Step 2: discretize triggering function and learn Quasi-conditional intensity. In this step, we choose piecewise constant function as the approximation to the true triggering function for the aggregated sequence. This means we will discretize the time horizon into small intervals (which we call bins) and estimate a "weight" on each interval. In practice, the time horizon we discrete is truncated on [0, T 0 ] and discretized into finitely many bins, since it is unnecessary to estimate infinite number of weights on infinite time horizon. More specifically, we estimate g 0 from the following class: G Here, 0 = δt 0 < δt 1 < · · · < δt n0 = T 0 and each bin B k = (δt k−1 , δt k ] has length ∆t k = δt k − δt k−1 (k = 1, 2, . . . , n 0 ). We apply Probability Weighted Histogram Estimation Marsan and Lengline [2008], Fox et al. [2016] to learn the weights g k on each bin B k , triggering magnitude α and background intensity µ. Details of this estimation is deferred to Appendix B.
Remark. Our Quasi-conditional intensity form satisfies the model assumption in Fox et al. [2016], which guarantees the non-parametric stochastic declustering algorithm is an EM algorithm. It will maximize a lower bound on the Quasi-log-likelihood function, which is in fact the complete-data Quasi-log-likelihood function derived by Veen and Schoenberg [2008]. Thus, it outputs MLE of Quasi-log-likelihood function (Quasi-MLE). See Appendix B for more details.
Step 3: compute GS statistic. Here, we call the singleton that we want to test a sub-model. We call the Quasi-parameter space under H 0 sub-model Quasi-parameter space and denote it by Θ 0 . Similarly, Θ is the full model Quasi-parameter space, or rather, Quasi-parameter space under H 1 . Under our proposed approximation class (2), the Quasi-conditional intensity has a parameterization θ = (µ, α (1) , g 1 , . . . , g n0 , α (2) , g 1 , . . . , g n0 ), where µ µ (1) + µ (2) . The full model Quasiparameter space is given by where d = 3 + 2n 0 . Note that the second constraint makes sure the Hawkes process is stationary and ergodic. Let n0 ) be the Quasi-parameter of the triggering function of Hawkes process z (z = 1, 2). The sub-model Quasi-parameter space is Denote the number of constraints (we'll see later it's in fact degree-of-freedom of our test statistic) r = 1 + n 0 = dim Θ − dim Θ 0 , the null hypothesis H 0 : θ 0 ∈ Θ 0 can be re-expressed as n0 ) . Note that here θ 0 is the parameter that maximizes the Quasi-likelhood and makes Quasi-conditional intensity a good approximation to the unknown true one. We consider a test: and the following test statistic: Definition 1 (GS statistic). Suppose the past sample trajectory is H T . Denote where H exists and has full row rank r. Then, the Generalized Score (GS) test statistic is given by where θ QM LE ∈ Θ 0 is Quasi-MLE under null hypothesis and Σ −1 is given by: Remark. Later, we will show T Σ −1 is a consistent estimator of inverse of covariance matrix of S T ( θ QM LE )/ √ T . Moreover, we give a closed-form expression for GS T in Appendix C.
Based on our testing procedure for two single data sequences above (steps 1 ∼ 3), we state a more general version for two sets of data sequences. Note that our theoretical results in the next section can be easily generalized to this general scenario, and we state them under two single data sequences scenario to ease the burden of notations.
Step I Mix D 1,i and D 2,i to get the aggregated sequence D agg i (i = 1, . . . , L).
Step II Apply Probability Weighted Histogram Estimation to learn Quasi-MLE.
Step III Repeat the following procedure for K times: randomly shuffle the order of sequences in the D 1 and repeat step I to get a different set of aggregated sequences, from which we randomly choose N sequences to calculate one GS T .

Theoretical Analysis
We present the asymptotic performance of our GS statistics by establishing a novel connection with classic results in statistics for MLE under model mismatch (Quasi-MLE in White [1982]). The asymptotic properties of MLE for Hawkes process without model mismatch are well studied in Ogata et al. [1978] and can be applied to a testing procedure like likelihood ratio test, score test and Wald test to get their asymptotic behaviors. In this section, we prove a generalization to model mismatch case for Hawkes process. The proofs of all theoretical results are in Appendix D.
Lemma 1 (Asymptotic properties of Qusai-MLE). Let θ QM LE and θ QM LE be Quasi-MLE under H 0 and H 1 , respectively. For piecewise constant triggering function family (2), Quasi-MLE satisfies the following asymptotic properties: (i) Convergence to θ 0 almost surely. That is, when T → ∞, (ii) Asymptotically normality. Define then, when T → ∞, we will have: Under H 1 : (iii) We also have asymptotically normality of the Quasi-score function, no matter under H 0 or H 1 : The asymptotic covariance matrix is inverse of the Fisher Information Matrix (FIM) I −1 (θ * ) for MLE and FIM I(θ * ) for score function in Ogata et al. [1978]. Moreover, by Theorem 1 in it, one can verify Information Matrix Equivalence Theorem in White [1982] still holds for stationary point process, i.e. θ 0 = θ * and A(θ 0 ) = B(θ 0 ) = I(θ 0 ) hold if and only if the model is correctly specified. Thus, it is easy to see our results simplify to the form in Ogata et al. [1978] in the absence of model mismatch. Though the asymptotic covariance matrix of Quasi-MLE is no longer inverse of the FIM, it can still be estimated consistently. Besides, our results can also be viewed as an extension of White [1982] to a stationary point process case.
Theorem 1 (Asymptotic null distribution of GS T ). Under the null hypothesis H 0 : θ 0 ∈ Θ 0 ⊆ Θ, the Generalized Score (GS) test statistic has an asymptotic χ 2 distribution. More specifically, Remark. Note that here the degree of freedom is r = 1 + n 0 ≥ 2 and depends on the discretization. For example, if one discretizes the time horizon [0, T 0 ] into 10 bins, then r = 1 + 10 = 11.
The asymptotic power function in shown is Figure 2. For more details of Marcum-Q-function, one can refer to Appendix E.

Numerical experiments
In this section, we present numerical results of our GS statistics to (1) validate the asymptotic property of our method by three simulation experiments; (2) demonstrate goodness-of-fit test for synthetic and real data.

Asymptotic property validation
To validate Theorems 1 and 2 presented in Section 4, we conduct three simulation experiments on a synthetic data set. To obtain reliable results, we repeat our experiments on five sub-data sets generated from Hawkes process defined in (1) with 1,000 sequences for each, where µ = 20 and an exponential triggering function φ(t − t i ) = αe −10(t−ti) , t i < t is adopted. The magnitude factor α in each sub-data set is from {1.5, 2, 2.5, 3, 3.5}, respectively.
The Q-Q plot in Figure 3 (a) shows that the GS statistic follows the χ 2 distribution, which confirms Theorem 1; Figure 3 (b) visualize the mean (red line) and the error bar (green bars) of each testing point for the GS statistics over different sample size N . As we can see, the GS statistics tend to be linear in sample size under H 1 , which matches the same theoretical results shown in our power study in Theorem 2 and validate that our asymptotic distribution analysis is reasonably accurate; Figure 3 (c) is the ROC Curve, where the true positive means we correctly accept the null hypothesis H 0 (identify two sequences are from the same distribution). The false-positive means we mistakenly accept the null hypothesis H 0 (identify two sequences are from the same distribution which they are not). We plot the same ROC plots with a different number of samples N used in calculating the GS statistics. We observe that the GS statistics show ideal classification performance when N = 100 (AUC is approximately 1). The testing procedure basically follows Section 3.2. We choose K to be 20, 5, 150 for three experiments, respectively. Details are deferred to Appendix F. In short, based on the following results, we have confirmed (a) the χ 2 null distribution; (b) the score is linear in sample size under H 1 ; (c) the consistency of the proposed test. In general, our proposed test can differentiate the subtle difference in triggering function accurately. We also conduct similar experiments for power triggering functions to validate our method is model free. The results are similar and we defer them to Figure 6 in Appendix F due to space limitation.

Goodness-of-fit test
We perform the complete testing procedure using our proposed GS statistics for four commonly used time series models on various synthetic and real data sets to demonstrate the testing power of our goodness-of-fit approach.
The experimental setting is described as follows. For synthetic experiments, we generate 5, 000 sequences for each data sets, which are generated from the Hawkes process (µ = 10) defined in (1) with different types of triggering functions: (a) exponential (Exp):  We follow the exact testing procedure described in Section 3.2 with N = 200, K = 5; we choose n 0 = 15 for Exp and Matern data and n 0 = 13 for MIMIC and MEME data. We report the mean of scores and the likelihood of fitting the model in Table 1. We conclude that our proposed goodness-offit test can differentiate models under different settings. In particular, the GS statistics can be used as a ranking criterion. More specifically, the parametric models Exp GD and Hist EM achieve lower scores (better performance) on synthetic data sets comparing to NHP and LSTM, since the parametric assumptions of the parametric models (e.g., the additivity in triggering effects) are consistent with the Hawkes process used in generating synthetic data. In the contrast, NHP performs better on real data sets, including MIMIC and MEME, where dynamics between events are more complex and difficult to be captured using parametric models. To make the numerical results more credible, we also present the corresponding likelihood in Table 1 as an indicator showing how well the data are fitted by the model (higher likelihood the better data is fitted). It has been shown that the likelihood result agrees with our GS statistics. Moreover, we also show that as a deterministic time series model, LSTM is difficult to compete with other baselines.
Goodness-of-fit for 911 call data. To demonstrate the use of our test statistic as a diagnosis tool for the goodness-of-fit of generative models, we test on 911 call data in 2017 provided by the Atlanta Police. The Atlanta Police Department divides its operation region into 78 beats, so we use this to partition the spatial region and consider a non-homogeneous point process generates sequences in each beat. We first consider police events data in each beats in one day as a sequence, and for each beat fit generative model using NHP and Exp GD. Then we calculated the value of the test statistic for each beat. The experiment configurations are as follows: N = 20, K = 1, n 0 = 12. The results are presented in Figure 4. Clearly, the generative model has different goodness-of-fit in each beat. Also, the two generative models have different patterns in their goodness-of-fit over space. Note that we do not know the ground truth. This example demonstrates that our tools provide a convenient and flexible diagnosis tool for the goodness-of-fit for generative models in practice. Each polygon in the map represents a police beat in Atlanta. The color depth represents the level of the test score. Lighter color means a smaller discrepancy between the generated data and the real data. Overall, NHP outperforms Exp GD on this data set in most of beats.

Broader Impact
With discrete event data, notably self-exciting discrete event data, being ubiquitous nowadays, modelling of such data to promote understanding of the inherent dynamics among those events is drawing more and more attention, especially in the fields like seismology, criminology which directly affect human life and property. Given that numerous sophisticated machine learning algorithms for capturing the interactions among those events have been developed so far, knowing how to rank those models and select the best one should have received more attention than developing a new one. In particular, choosing the better model for the triggering function for Hawkes process is essential for improvement in forecasting algorithms based on it, which in turn will benefit the society (e.g. guiding earthquake evacuation). For example, it has been found that spatial power law triggering function is much more effective than normal kernel and the former one will improve the earthquake forecasting algorithms Reinhart et al. [2018]. To the best of our knowledge, our proposed test is the first to differentiate the subtle difference in triggering function and can be of great help in selecting the best model in real study. We should note that this testing procedure is limited to Hawkes process, i.e. linear conditional intensity and additive self-exciting (or correcting) effect as specified in (1). Nevertheless, this is already general and powerful enough to capture the desired interactions among those events and help decision-makers to make those life-concerning decisions.
Erik Lewis, George Mohler, P Jeffrey Brantingham, and Andrea L Bertozzi. Self-exciting point process models of civilian deaths in iraq. Security Journal, 25 (3)

A Mixing of two Hawkes processes
We first present a useful lemma, which provides the proof for full model case (i.e. under H 1 ). Another equivalent definition of conditional intensity λ(t|H t ) for a counting process {N (t) : t ≥ 0} with associated history H t (t ≥ 0) is We will make use of this definition to prove the following lemma.
Lemma A.1. Suppose we have n Hawkes processes {N z (t) : t ≥ 0} (z = 1, 2, . . . , n) with conditional intensity specified by (1). Define the mixing to be N (t) = n z=1 N z (t). The conditional intensity of mixing of n Hawkes processes is sum of those n conditional intensities. That is, The occurrence of next event may be triggered by either history of Hawkes Process 1 or Hawkes Process 2 and thus it is dependent on the joint history.

Events in Hawkes Process 1
Events in Hawkes Process 2 Figure 5: Illustration of mixing of two Hawkes processes N (t). Given the past sample trajectory, the upcoming event of N (t) may (1) come from background poisson process of Hawkes process 1 or 2 OR (2) be a offspring of history H 1,t or H 2,t . The grey dashed line in the figure illustrated scenario (2).
Proof of Lemma A.1. We prove by the definition of conditional intensity. For any non-negative integer m ∈ Z + , denote :: m = (m 1 , . . . , m n ) and M = { :: m | m 1 + · · · + m n = m, m i ∈ Z + }, Case 1: When m > 1, it is easy to see P(N (t + h) − N (t) = m|H t ) = o(h), since either there are at least two m i 's ≥ 1 or at least one m i ≥ 2.
Case 2: When m = 1, there will be one and only one of all m i 's taking value 1 and the rest will be all zeros. Thus, we have Case 3: When m = 0, all m i 's will be zeros and we will have Let λ(t|H t ) = n i=1 λ i (t|H i,t ), and we will find out this is the conditional intensity for N (t).
Proof of Proposition 1. We can see under the alternative hypothesis, the result directly follows Lemma A.1. Under null hypothesis, by Lemma A.1, it is easy to show N (t) defined in Proposition 1 has intensity By the definition of Hawkes Process in Section 2, we can see the mixing of two Hawkes processes under H 0 is still a Hawkes process. We complete the proof.

B A non-parametric estimation of the Quasi-conditional intensity B.1 Probability Weighted Histogram Estimation under null hypothesis
Define the branching structure as follows: probability event i is triggered by event j, i > j probability event i comes from background, i = j 0, i < j Apparently, we want to estimate the Quasi-background intensity from background process only and Quasi-triggering function from the triggered events only. Instead of using a hard-threshold indicator, Zhuang et al. [2002] used a stochastic declustering procedure to separate the background events from triggered ones by assigning each event a weight, or rather the probability that this event comes from background or is direct offspring from an individual ancestor. Then, we can use a probability weighted estimator to estimate Quasi-background intensity and Quasi-triggering function. The algorithm is as follows: Assume we have estimated branching structure p (v) ii at iteration v, then we can estimate the Quasibackground intensity as follows: For the Quasi-triggering component, as we assume g to be a p.d.f., we can estimate the magnitude of triggering effect from triggered events only: For the temporal component in the Quasi-triggering function, for each bin (as we discretize in (2)), we estimate its parameter from those triggered events which falls into that bin, i.e.
After estimating the Quasi-conditional intensity function, we update the branching structure. More specifically, for i > j: We summarize the algorithm as follows: Algorithm 2 Probability Weighted Histogram Estimation of Quasi-log-likelihood under H 0 Initialize: choose stopping critical value (e.g. 10 −3 ), initialize p ij and set p 1. Estimate Quasi-background rate µ as in (3). 2. Estimate Quasi-triggering components magnitude α and temporal g(t) as in (4) and (5). 3. Update probabilities p (v+1) ij 's as in (6) and (7). 4. v ← v + 1 end while B.2 Derivation of Probability Weighted Histogram Estimation as EM-type algorithm under null hypothesis Fox et al. [2016] demonstrated that algorithm 2 is an EM-type algorithm under (2) by using complete data log-likelihood. However, they did not explicitly show the E-step also maximizes the complete data log-likelihood (or we should say complete data Quasi-log-likelihood here since Fox et al. [2016] assumed the ground truth takes piecewise constant form (2)).
At high level, we will first lower bound the Quasi-log-likelihood and then show that the algorithm iterates between maximizing this lower bound w.r.t. branching structure (p ij 's) and w.r.t. the Quasi-conditional intensity (Quasi-background rate µ, Quasi-triggering magnitude α and temporal Quasi-triggering function g). First recall the Quasi-log-likelihood function under H 0 : We can simplify the last term above by using integral approximation of Schoenberg [2013]: Thus we can ignore the last term when maximizing the log-likelihood function.
Next, we lower bound the first term in the Quasi-log-likelihood function by Jensen's inequality: where p ij 's satisfy i≥j p ij = 1. Then we can get a lower bound on the approximation of Quasilog-likelihood under the piecewise constant parameterization: Denote this lower bound by (θ). We maximize this lower bound under the following constraints: n0 k=1 g k ∆t k = 1, (g(t) is a p.d.f.) i≥j p ij = 1, (p ij 's are probability weights) By adding Lagrange multipliers, this is equivalent to maximizing the following objective: M-step: By taking first order derivative w.r.t. µ and setting it to zero, we will have: Solving for µ and we will get which is the same with the update in step 1 in Algorithm 2. This means when we have p ij 's at iteration v, the update in step 1 in Algorithm 2 leads to a larger Quasi-log-likelihood value. Similarly taking derivative w.r.t. α and setting it to zero leads to the update in step 2: Solving for g k and c 1 by some simple algebra and we will get the update for g k when we have p ij 's at iteration v as follows: Repeat the similar procedure, we will get: By the first two equations we have Plug this back into the last equation and we will get the update in step 3 in Algorithm 2. Thus, we validate Algorithm 2 as an EM-type algorithm.

C Explicit form of GS statistic
Note that g (z) (u) = n0 k=1 g (z) k 1 B k (u). To simplify the explicit expressions, we first define the following notations: Here, G(i, z ; z) represents the triggering effect of events in process z to i−th event in process z .
First order partial derivatives: Here, we get the explicit expression for S T (θ) and A T (θ).
Second order partial derivatives: We first show that the assumptions in Ogata et al. [1978] hold for our Quasi-conditional intensity function.
(B) The Quasi-conditional intensity function we consider here is actually linear w.r.t. the parameters, then it is arbitrarily order continuous differentiable (i.e. smooth) and bounded within any compact set in the Quasi-parameter space. Assumptions (B1) ∼ (B7) hold trivially.
(C) By (2), the Quasi-temporal triggering function is truncated on [0, T 0 ], which means and complete data conditional intensity function λ(t|H −∞,t ) will be exactly the same with λ(t|H 0,t ) as long as t > T 0 . Since Assumptions (C1) ∼ (C4) only require stochastic approximations of λ(t|H 0,t ) to λ(t|H −∞,t ) when t goes to infinity, it is easy to see those assumptions are satisfied.
Next, since our parametric form (2) is only approximation to the true one, we need to slightly modify the theoretical results in Ogata et al. [1978] for our Quasi-MLE. Here we will not mention theorems or lemmas that we do not need to modify under model mismatch (except that we should keep in mind that the "true" parameter in Ogata et al. [1978] is understood as the maximizer of Quasi-likelihood) and it is easy to verify those theoretical results (from the beginning to Theorem 5) by just following the proof therein.
Before we proceed to the proof, we should note that the Quasi-MLE is θ QM LE under H 0 and θ QM LE under H 1 . Under H 1 : θ 0 ∈ Θ 0 , the estimator θ QM LE is obtained using the full model conditional intensity 1 instead of 0 . The estimation is given in Algorithm 3 in Appendix F.
For simplicity, we denoteθ QM LE to be θ QM LE and θ QM LE under H 0 and H 1 , respectively. That is, Modifications on Theorem 1. Here θ 0 is not the true parameter of the true conditional intensity function. Instead, it is the maximizer of Quasi-log-likelihood, i.e. our approximation to the true log-likelihood function. By the definition of θ 0 and stationarity of the process, the first result still in this theorem still holds: However, the second result does not hold unless our approximation is indeed a correct specification of the model. More specifically, in general, where λ * is the correct parametric form and typically unknown in practice.
Thus, we have Using our notation, this can be re-expressed as A(θ 0 ) = B(θ 0 ).
Modifications on Theorem 2. The convergence in our case is much stronger. By following the proof in Fox et al. [2016], the convergence in probability comes from Assumptions (C), where the convergence in the stochastic approximation is only in probability sense. However, we just show that the stochastic approximation holds for every sample path as long as t > T 0 based on our parameterization (2) that the Quasi-temporal triggering function is truncated, i.e. our convergence is in almost surely sense. Thus, we have: Modifications on Theorem 4. Since A(θ 0 ) = B(θ 0 ), the convergence result should be This is because where the first equality comes from definition and stationarity of the process.
Modifications on Theorem 5. By the proof of this theorem one can reach this result: Again, since A(θ 0 ) = B(θ 0 ), the asymptotic covariance matrix is not B −1 (θ 0 ) and that's the modification here. Besides, the asymptotic χ 2 distribution of log-likelihood ratio does not hold because of the model mismatch.
Here, we complete the proof.

D.2 Proof of Theorem 1: asymptotic distribution under null hypothesis
This proof is highly involved. To help better understand this proof, we first provide a high level sketch on why our GS statistic follows a χ 2 distribution.
By adding Lagrange Multiplier ζ T , we can derive that θ QM LE satisfies: Following idea in Boos [1992], we can use Taylor expansion to expand S T (θ 0 ) about θ QM LE and h( θ QM LE ) about θ 0 (note that we have h(θ 0 ) = 0 under H 0 ): Note that by our notation ∇h(θ) = H(θ). Since h(θ) is linear in θ, its gradient is a constant matrix and we can denote H = ∇h(θ).
Pre-multiply the first equation above by H HB −1 The matrix in the middle of RHS is a projection matrix for the column space of B − 1 2 T ( θ QM LE )H , and from (8) is already in this space. This means the RHS is exactly S T ( θ QM LE ) and we will get: Rewrite GS statistic as By Lemma 1, one can verify S T ( θ QM LE )/ √ T has a asymptotic normal distribution with T Σ −1 being a consistent estimator of generalized inverse of its asymptotic covariance matrix.
Since H is of rank r, we verify that GS T ∼ χ 2 r .
Next, we present a more rigorous proof following the method in White [1982].
Proof. We first state some useful results: By the almost surely convergence of Quasi-MLE (modifications of Theorems 2 and 5 in Ogata et al. [1978]), we have that The modification of Theorem 1 in Ogata et al. [1978] can be re-expressed as S(θ 0 ) = 0.
The modification of Theorem 4 in Ogata et al. [1978] can be re-expressed as follows where S T is the Quasi-score function (i.e. first order gradient of Quasi-log-likelihood function).
Under null hypothesis, the asymptotic χ 2 distribution of GS statistic under model mismatch (e.g. Theorem 3.5. in White [1982] and Section 4.2. in Boos [1992]) can be extended to temporal Hawkes process.
Validation of asymptotic properties in Section 5.1: (a) For each α, generate L data sequences as D 1 and another L data sequences as D 2 ; (b) Generate L data sequences from α = 1 as D 1 and another L data sequences as D 2 from α = 4; (c) Use the first pair of data sequence set in (a) (corresponding to α = 1.5) as positive sample and data sequence set in (b) as the negative sample.
The experiment configurations (initialization) are as follows: L = 1, 000, n 0 = 14 and end points for those bins are (0, .04, .08, .12, .16, .2, .26, .32, .38, .45, .55, .65, .75 Goodness-of-fit in Section 5.2: D 1 is chosen to be the testing data and D 2 is generated from the model fitted on the training data. For 911 call data, L = 364 and we use the first 200 sequences to as the training data to fit the model and the rest 164 sequences as D 1 . Then we generate 164 data sequences as D 2 to perform the testing procedure. We choose N = 20, K = 1 and use (0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.5, 1) as end points for bins.

F.2 Additional experiments
To show our proposed non-parametric method is able to detect subtle difference in triggering functions under different settings, we use different synthetic data to validate our theoretical results.

G Probability weighted histogram estimation under alternative hypothesis
Under H 1 , the triggering mechanism is more complex compared to univariate Hawkes Process, since each event can be either from the background, direct offspring from an individual ancestor in Hawkes Process 1 or Hawkes Process 2 and the triggering effects of events in two different processes are different.
We denote the branching structure matrix P (z)(z ) ∈ R (N1+N2)×(N1+N2) (z, z ∈ {1, 2}). The element in i−th row and j−th column is defined to be the probability that event i in process z is triggered by event j in process z if either i = j or z = z (case 1) or the probability that event i in process z is a background event if i = j and z = z (case 2). That is, probability that event i in process z is triggered by event j in process z , case 1 probability that event i in process z is a background event, case 2 Note that the probability is zero when t (z) i ≤ t (z ) j , which means the event that happens earlier in the process cannot be triggered by those which happen later.
As discussed above, we focus on differentiating difference in triggering effect. Thus we estimate the sum of two background intensities from all background events: For the triggering components, we estimate the magnitude for process z (z = 1, 2) using events from aggregated data triggered by process z and estimate the temporal triggering density function from those events which fall into the corresponding bin. Note that we have This is because the probability will be zero if t (z ) i ≤ t (z) j as discussed above. Thus, for z = 1, 2 and k = 1, . . . , n 0 , the estimators can be expressed as Nz j=1 P (v) (z )(z) ij framework above, we first use integral approximation of Schoenberg [2013] to approximate the Quasi-log-likelihood function and then lower bound it using Jensen's inequality: where z = z. Note that the term −N 1 α 1 − N 2 α 2 comes from integral approximation. Add Lagrange multipliers and we will get the following objective function: Then, by taking first order derivatives and setting them to zero we can validate Algorithm 3 as an EM-type algorithm.