{'fraction_non_alphanumeric': 0.047596470433293826, 'fraction_numerical': 0.015487949427103912, 'mean_word_length': 4.0286092715231785, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A test for dependence between two point processes on the real line

Many scientific questions rely on determining whether two sequences of event times are associated. This article introduces a likelihood ratio test which can be parameterised in several ways to detect different forms of dependence. A common finite-sample distribution is derived, and shown to be asymptotically related to a weighted Kolmogorov-Smirnov test. Analysis leading to these results also motivates a more general tool for diagnosing dependence. The methodology is demonstrated on data generated on an email network, showing evidence of information flow using only timing information. Implementation code is available in the R package `mppa'.


Introduction
Testing for dependence between two point processes is a long-standing statistical problem. When the two processes are on the real line, usually representing time, the points are usually interpreted as event times (or, simply, events). Scientific questions then often revolve around detecting: (a) Triggering: the occurrence of an event in A temporarily increases the rate of events in B.
(b) Correlation: the rate of events in B is locally increased around events in A.
(c) Inhibition: the occurrence of an event in A temporarily decreases the rate of events in B.
(d) Anti-correlation: the rate of events in B is locally decreased around events in A.
Statistical methods to detect these effects have received decades of attention in the field of neurophysiology. The activity of a neuron is often recorded as a sequence of 'spike' times, called a neuronal spike train, which is often treated as a realisation of a point process on the real line. Comparing trains that are simultaneously generated by different neurons can shed light on how they are connected and, more generally, how information is processed in the nervous system. The literature in this field is relatively mature, for example a very highly cited paper by Perkel et al. (1967) proposed to test for interaction on the basis of histograms of the times from A to B events. A number of difficulties were recognised, namely, that the proposed tests are not exact even when A and B are independent homogeneous Poisson processes and that the histograms can be misleading if the processes have an oscillatory component of similar frequency (causing spurious correlation when the processes are 'phase-locked' by chance). Since then, a number of model-based approaches were developed for this problem, notably in a series of papers by Brillinger (Brillinger et al., 1976;Brillinger and Segundo, 1979;Brillinger, 1988bBrillinger, ,a, 1992. In a more general context Ripley (1976Ripley ( , 1977 introduced the so-called K-function to measure second-order dependence between point processes defined on a topological space. From this work a number of articles followed, typically motivated by ecological or biological applications, adapting ideas to two-or three-dimensional settings. Lotwick and Silverman (1982) proposed to use an estimate of K to test for interactions between different types of points on a plane. Berman (1986) presented some results for when B is a spatial point pattern and A is an arbitrary stochastic process. Back in a one-dimensional setting Doss (1989), also motivated by a neurophysiological application, proved that a one-dimensional version of the estimated (cross)-K function is asymptotically normal.
The problem has now resurfaced in the analysis of network data, for example traffic generated on a computer network, messages and other connections made on social networks, mobile communications, email networks, the web, collaboration networks (e.g. in academia, music or film) and more. Such data can often be represented as a graph with point processes (e.g. communication times) occurring on every edge (e.g. a pair of computers). In being able to measure dependence between the events generated by edges or nodes, there is great potential to better understand information flow, discover new correlations and develop better network models. Recent approaches include Blundell et al. (2012), where email reciprocation is modelled using the mutually exciting point process models developed by Hawkes (1971a,b), and Perry and Wolfe (2013), developing a framework for modelling point process networks based on a version of Cox's proportional intensity model. In this approach effects such as A triggering B are handled by including past events of A as covariates for B.
This article seems to be the first to propose an exact (generalised) likelihood ratio test for association between two one-dimensional point processes. We assume we are given a model for B, which captures its statistical behaviour under the null hypothesis, for example encapsulating any seasonality, changepoints or drift. Then we test for a multiplicative effect on the intensity of B within a certain interval following or surrounding every A event. The multiplier and interval width are fit by maximum likelihood estimation, leading to a (generalised) likelihood ratio test for the corresponding effect. Depending on whether the alternative hypothesis assumes an increase or decrease in the intensity of B events around A events, and whether the intervals are centered on the event times or occur directly afterwards, the same procedure can be used to test for each of the effects listed in (a-d).
The tests obtained are easily implemented; their p-values are available in an exact (albeit recursive) form and they are shown to be asymptotically related to a weighted Kolmogorov-Smirnov (K-S) test. A key insight is noticing a duality with finding a changepoint in a homogeneous Poisson process (Lemma 1), after which mathematical considerations are greatly simplified. A by-product of the lemma is a new diagnostic tool for analysing dependence between point processes.
The remainder of this article is organised as follows. In Section 2 we focus on testing for events in A triggering events in B under a non-homogeneous Poisson null hypothesis. We give the test and derive its finite sample p-value. Section 3 proposes a number of adjustments to improve power against some alternatives and extends the methodology to detecting other forms of dependence. In Section 4 analysis is extended to more general point processes on the real line. In Section 5 an asymptotic analysis finds a connection between our test and a weighted K-S test of a set of auxiliary variables constructed during the procedure. Finally in Section 6 the test is applied to detecting reciprocation and forwarding behaviour in an email network (the Enron email dataset).

A fixed duration excitation model
Let A and B be two simple point processes on the real line observed simultaneously from the first event time of A up to an observation end time L. Neither process is explosive, so that the observed event times of A form a finite set A = {a 1 < . . . < a m }, m ≥ 1, and the event times of B form a finite set B = {b 1 < . . . < b n } where b 1 ≥ a 1 and n ≥ 0. Without loss of generality, let a 1 = 0. We make no further assumptions on A, and treat its event times as given when modelling B.
In this section we consider the special case where B|A is a non-homogeneous Poisson process. Relaxations of these conditions are considered in Section 4. Given A, the process B is assumed to have a deterministic, bounded, Lebesgue-measurable intensity function λ B (t). Heuristically λ B (t)dt = EdB(t).
We propose the following fixed duration, proportional excitation model for B: where a(t) is the most recent event in A occurring at or before t, τ > 0, λ 1 ≥ λ 2 ≥ 0 are unknown parameters, and r is a known bounded non-negative Lebesgue-measurable function satisfying L 0 r(v)dv = 1. Model (1) leads to a test with a very straightforward interpretation: is the relative proportion of events within time τ of an event in A higher than can be explained by r alone? To formalise this idea, we shall test the following null hypothesis: A natural test statistic is the generalised likelihood ratio where is likelihood of B under model (1). Note that under H 0 , τ has no real importance and, accordingly, (B; τ, λ 1 , λ 2 ) is functionally independent of τ when λ 1 = λ 2 .
Remark 1. Probably the most common choice for r will be the constant 1/L, in which case B is a homogeneous Poisson process under H 0 . If the application makes this assumption unrealistic, a more informed choice of r will not only bring the null behaviour of the test closer to its nominal distribution, derived later, but may lead to a gain in power under the alternative, for example if B responds to A despite being relatively inactive under H 0 . If r is unknown our (heuristic) recommendation is to use an estimate (as we do in Section 6).
Remark 2. If B was observed for some time before a 1 , a period that is ignored above, it could be modelled as being generated by λ 2 r(t) for t < a 1 . Doing so only requires some simple modifications to the results below.

Computing the test
Let ρ(X) = X r(v) dv, where X is any set that can be constructed by a countable union of intervals. ρ can be seen as an artificial measurement of time intervals that compensates for B's activity levels under H 0 . In fact, ρ is a Borel probability measure on [0, L). The likelihood of Model (1) can be written (Daley and Vere-Jones, 2003, p. 232) where T(τ ) = {t : t − a(t) ≤ τ } is the union of all intervals during which B was intensified due to the events of A and K(τ ) = #{b i ∈ T(τ )} is the number of events that occurred within them. The likelihood is maximised as follows. Let u 1 ≤ ... ≤ u n be the order statistics of ρ{T(b i − a(b i ))}, i = 1, . . . , n. u k can be interpreted as the effective proportion of time that B was intensified if τ is equal to the kth smallest response time. The maximum likelihood parametersτ ,λ 1 ,λ 2 are found within the following procedure: Procedure 1 (Computation of T given u 1 , . . . , u n ). For k = 1, . . . , n, let The maximum of (3) is a monotonic function of k . This can be shown by a straightforward argument, given in Appendix A. Note that k = n satisfies u k ≤ k/n so thatk is always defined.
Since the number of events of B can be equally well explained under the null as under the alternative, it is natural to condition on the value of n. The denominator of (2) is a multiplicative term that depends only on n and can then be ignored. Hence any monotonic function of the numerator is admissible as a test statistic, and we use T = k .

A simple reformulation
The following lemma establishes a duality between testing for model (1) and detecting a Poisson process change-point, and is the key observation of this article.
This result, although it is very simple, is new to the best of the authors' knowledge. However similar results exist, notably Theorem 1 in Loader (1991) on finding a changepoint in the hazard rate of n independently and identically distributed (IID) survival times. Lemma 1 and its proof provide a number of fascinating insights into the testing problem.
Testing for a Poisson process change-point. Revisiting Procedure 1, we can see that T only depends on the data through {u 1 , u 2 , . . . }. More than that, T is the generalized likelihood ratio test for model (4) against a homogeneous Poisson null hypothesis, A diagnosic tool. Conditional on n, u 1 , . . . , u n are ordered Uniform random variables under the null hypothesis, whereas under the alternative they should be, roughly speaking, more concentrated towards 0. To assess dependence the u i can be inspected without having model (1) in mind, for example using goodness-of-fit tests against uniformity or visual diagnosis tools (e.g. a plot of the empirical cumulative distribution function (CDF) of u i compared to y = x). In fact, in Section 5 we show that T is asymptotically related to a weighted Kolmogorov-Smirnov test.
'Chopping and stretching'. The proof essentially relies on superposing segments of B and obtaining the intensity function of the superposed process by adding the intensities contributed by each of the relevant segments. This is possible because of the Poisson assumption. Then time is 'stretched' according to that intensity, leaving only the effects of the unknowns, λ 1 , λ 2 and τ . This intuition invites us to consider other ways of superposing B based on A, and motivates a number of tests to follow.
Proof. Notice that K(y) = #{b i ∈ T(y)} defines a point process over [0, max{a i+1 − a i : i = 1, . . . , m}), temporarily defining a m+1 = L. Its increments are obtained from the event times of B by superposing the segments [a i , a i+1 ), aligning to the left. The intensity of K(y) is therefore h(y) = {t:t−a(t)=y} λ B (t) and its compensator is If F is a continuous non-decreasing function on a sub-interval X of R with image Y , we define F −1 (y) = inf{x ∈ X : F (x) ≥ y} for y ∈ Y . Because ρ is continuous, ρ{T(y)} is a continuous non-decreasing function µ, say, of y. Then µ −1 is right-continuous and nondecreasing with jumps at an at most countable set of values of x corresponding to intervals where µ is constant (Daley and Vere-Jones, 2007, p.420) is the qth smallest response time (and the qth event time of K). Using the well-known time change theorem (Daley and Vere-Jones, 2007, p.421 is homogeneous Poisson with intensity λ 1 over that range. By a similar argument U (x) has intensity λ 2 over (ρ{T(τ )}, 1).

Distribution of T
Let F n (x) = P(T ≤ x|n) denote the distribution function of T under the null hypothesis, conditional on n. To simulate T from F n we can generate n uniform variables independently, set their sorted values to u 1 ≤ · · · ≤ u n , and compute T through Procedure 1. From this the moments and critical values of T can be approximated by simulation.
On the other hand the p-value for an observed test statistic t, P[T ≥ t|n], can be computed exactly. To this end notice that where o i is the solution for x ∈ (0, i/n] to which is obtained numerically. Plainly speaking o i is the smallest u i can have been if we know T = t. Worsley (1988) has made a similar observation in the context of testing for a changepoint in the hazard rate of IID survival times. (Although the test considered there is two-sided.) Various recursive formulas exist for computing the joint survival probability of n ordered Uniform variables, although many are unsuitable for computation because they involve differences of very large numbers (so-called catastrophic cancellation). A 'safe' option is the O(n 2 ) formula in Noé and Vandewiele (1968), as corrected in Noé (1972). This recursion is implemented in the R package corresponding to this article, 'mppa'. For the reader's convenience, it is also available in pseudo-code given below.
Algorithm 1 (The joint survival probability of ordered Uniform variables).

Adjustments and other types of association
In this section the test proposed in Section 2 is extended in two senses. First, we propose some adjustments to improve power against some more specific triggering alternatives. Second, conceptually similar methodology is developed to detect other sorts of dependence.

Time-limited τ
We may be interested in limiting τ to a maximum range, τ max say. This is useful if τ is expected to be small under H 1 , in which case restrictingτ ≤ τ max avoids 'wasting' power on testing for longer-term dependence.
To obtain a p-value of this test report p = 1 if T = 1. Otherwise compute o 1 , . . . , o n as in (6), and calculate using Algorithm 1. Note that the p-value of this test has a distribution with a discrete component: it is exactly one whenever no events occur within τ max . Nevertheless for any α ∈ [0, 1] we have P[p ≤ α] ≤ α under the null hypothesis, resulting in a conservative test.

Adaptive range
Consider the slightly modified alternative hypothesis In this model how long B is intensified following events of A is adapted to the null intensity of B: B has longer to respond if relatively inactive according to r.
This observation provides a testing procedure for model (7). A test statistic T is computed as in Section 2.1, using as inputs A * and B * instead of A and B, and using the Lebesgue measure for ρ, that is, ρ(X) = X 1 dx. Like the original test, T has distribution F n under H 0 conditional on n.

Other types of association
The following tests simply replace the input to Procedure 1, previously u 1 , . . . , u n , by a sequence that also satisfies (4). The output is a likelihood ratio test for a different form of association, but which has the same null distribution, F n , conditional on n.
Correlation. Relax the constraint a 1 = 0. Letã(t) be the closest event to t in A, which can now occur before or after t. In (1), replacing t − a(t) ≤ τ by |t −ã(t)| ≤ τ results in a model where the rate of events of B is locally increased around events of A. A likelihood ratio test is obtained as follows. Let C(τ ) = {t : |t −ã(t)| ≤ τ } and let v 1 ≤ ... ≤ v n be the order statistics of ρ{C(a i −ã(b i ))}, for i = 1, . . . , n. Compute T by inputting v 1 , . . . , v n to Procedure 1. By a similar argument to the proof of Lemma 1, we find that v 1 , v 2 , . . . are the event times of a point process following model (4). Thus T has distribution F n under H 0 , conditional on n.
Inhibition. In model (1), constraining λ 1 ≤ λ 2 and replacing t − a(t) ≤ τ by t − a(t) < τ (to simplify analysis), we obtain a model for inhibition. By calculations analogous to those of Section 2 we find that the likelihood ratio test T for this model can be computed by simply inputting x i = 1 − u n−i+1 , i = 1, ..., n in Procedure 1. As before T has distribution F n under H 0 , conditional on n.
Anti-correlation. Putting λ 1 ≤ λ 2 and replacing t − a(t) ≤ τ by |t −ã(t)| < τ we have a model where the rate of events of B are decreased around events of A. To test for this, use y i = 1 − v n−i+1 , i = 1, ..., n in Procedure 1.
Role reversal. For each of these forms of association analysis could conceivably be turned around by testing A given B. This only needs reversing time and swapping A and B in Section 2 (but note that an estimate of the null intensity of A is required). These role and time reversals result in different test statistics from those obtained by testing B, even though they target a similar type of association. It is an interesting question what circumstances could lead us to choose this approach.

Extension to more general point processes
In this section we consider relaxations of the non-homogeneous Poisson assumption in Section 2. Lemma 1 does not easily translate to a general point process if r in (1) is simply replaced by a conditional intensity (formally defined in Section 4.2). This also invalidates the results that follow, including the null distribution of the test. The difficulty is illustrated in a simple example in Appendix B, where B is a point process with just one point uniformly distributed on [0, T ), and A = {0, T /2}.

Independence conditional on n
A simple generalisation is possible if we can replace the intensity by a density over the observation region. Suppose that B can be generated by drawing n from some distribution and then placing the n event times independently according to some probability measure over [0, L). The non-homogeneous Poisson process is an example, as is the single-point point process mentioned above. Further remarks on the construction and properties of such processes are given in Daley and Vere-Jones (2003, Chapter 5). Analogously to model (1), suppose that the measure has a density where r is as before a bounded non-negative Lebesgue measurable function satisfying L 0 r(v)dv = 1, thereby defining ρ. If n has a Poisson distribution this model reduces to (1). The hypothesis test H 0 : λ 1 = λ 2 versus H 1 : λ 1 > λ 2 can be evaluated through a similar generalised likelihood ratio test: we first compute u 1 , . . . , u n as the order statistics of ρ{T(b i − a(b i ))}, i = 1, . . . , n and then T using Procedure 1. By straightforward modifications to the proof of Lemma 1, we find: Lemma 2. Under model (8), conditional on n, {u i } is a set of IID random variables with support on [0, 1) and density From this we establish that T given n also has null distribution F n under model (8).

Random time change
A much wider class of point processes can be conceived by allowing the intensity at time t, previously a deterministic quantity, to be dependent on the past. Formally, information available at time t is captured by an increasing family of σ-algebras F = {F s : 0 ≤ s < ∞} with F s ⊆ F t for 0 ≤ s ≤ t < ∞. If B is adapted to a right-continuous history F, by the Doob-Meyer decomposition there exists a right-continuous predictable process Λ(t) such that for every i ∈ N the stopped process This process is called the compensator for B. Under mild regularity conditions the conditional intensity λ(t) of B can then be defined via Λ(t) = t 0 λ(x)dx (Daley and Vere-Jones, 2007, p.358, p.367, p.390).
Note the following result, more general than that needed in the proof of Lemma 1, known as the random time change theorem (Daley and Vere-Jones, 2007, p.421): if B is non-terminating (there are infinite events as t → ∞) and Λ 0 is continuous then the process is a homogeneous Poisson process with unit rate. (Recall our definition F −1 (y) = inf{x : F (x) ≥ y} for a non-decreasing F .) Suppose then that under the null hypothesis B has a conditional intensity λ 0 . The random time change theorem allows construction of a processB : t ∈ [0, Λ 0 (L)), which, under the null hypothesis, is a stopped unit-rate homogeneous Poisson process. To compare like for like, letÃ(t) = A(Λ −1 0 (t)). Since we cannot get exact distributional results when we replace r in (1) by λ 0 a natural alternative might be to test forÃ triggeringB via the model whereã(t) is the most recent event inÃ occurring at or before t. A testT for this model is obtained by computing the corresponding u i and then implementing Procedure 1 replacing b i by Λ 0 (b i ), a j by Λ 0 (a j ), L by Λ 0 (L) and finally r(t) by the constant 1/Λ 0 (L).
Although conceptually similar to the test described in Section 3.2, where τ was adapted the null intensity of B, the dependence of the stopping-time Λ 0 (L) on B under H 0 makes inference more complicated. In particular Lemma 1 and the p-value computed in Section 2.3 no longer hold exactly, but may be sufficiently close approximations for practical use. Progress is possible if we allow L to be random, but this seems contrived.

Asymptotic distribution
In this section we demonstrate an asymptotic connection between our test and a weighted Kolmogorov-Smirnov (K-S) test. An especially surprising aspect of this result is that the weight function needed is exactly that which is most often used for K-S tests, e.g. in Chicheportiche and Bouchaud (2012). Note that classical asymptotic results about the likelihood ratio test do not hold here because of the discontinuity of the likelihood function at τ under H 1 (Müller and Wang, 1994, p.227).
Let y 1 , . . . , y n be independent replicates of an absolutely continuous random variable Y with support Y ⊆ R and CDF F , and consider their empirical CDF A generalized K-S test (Anderson and Darling, 1952) can be used to test whetherF is consistent with a hypothesized distribution F , for some weight function φ(x) ≥ 0, x ∈ [0, 1]. We need a one-sided, interval-restricted version of the above, for 0 < γ 1 < γ 2 < 1, defining G + [γ 1 ,γ 2 ] = 0 if the set above is empty. In the theory of K-S tests, it is often a consequence ofF being right-continuous and φ satisfying a specific property given in Anderson and Darling (1952, Eq 2.8) that attention can be restricted to the value of F at the observed points. This is not the case here: the presence of γ 1 and γ 2 forces us to define G + [γ 1 ,γ 2 ] as (11). Analogously, let T [γ 1 ,γ 2 ] = max{ i : u i ≤ i/n; γ 1 ≤ u i ≤ γ 2 }, now setting T [γ 1 ,γ 2 ] = 1 if the set above is empty.
Let n, i → ∞ with i/n = p. Then if p ∈ [γ 1 , γ 2 ], u i is asymptotically normal with mean p and variance proportional to 1/n and therefore u i = p + O p (1/ √ n). Viewed as a function of u i , i has first and second derivatives 0 and p −1 (1 − p) −1 . Therefore by Taylor expansion we find i = 1 + (p − u i ) 2 /[2p(1 − p)] + o p (1/n) (Davis and Brockwell, 1991, Prop. 6.1.5). Hence Let S = {i : u i ≤ i/n; γ 1 ≤ u i ≤ γ 2 }. Heuristically ignoring the influence of o p (1) terms (of which there are a growing number with n), assume defining the RHS to be zero if S is empty. The absolute value was removed because i ∈ S guarantees u i ≤ p. Now, use p =F (y i ) and F (y i ) = d u i in the numerator andF (y i ) → p in the denominator to find that the RHS above is also the limiting distribution of G + [γ 1 ,γ 2 ] . To prove this more rigorously we would need a better understanding of the joint behaviour of the o p (1) terms (if arbitrary random variables x 1 , ..., x m are each individually o p (1) then max(x i ) is not o p (1) in general if m increases with n), perhaps reasoning along the lines of Donsker (1952) (on replacing a sequence or order statistics of a Uniform random variable by its limiting process in the asymptotic theory of K-S tests).
By simulation it seems that result does not hold in the limit γ 1 → 0, γ 2 → 1. This is not necessarily surprising -for instance it is noted in Chicheportiche and Bouchaud (2012) that with γ 1 = 1/(n + 1) and γ 2 = n/(n + 1) the asymptotic distribution of the two-sided version of G [γ 1 ,γ 2 ] still depends on n, and in fact the asymptotic theory of weighted K-S tests generally relies on φ being bounded over the unit interval (Anderson and Darling, 1952, p.196), which [p(1 − p)] −1 is not.

Example: information flow in the Enron email corpus
The Enron email corpus is a dataset that comprises emails sent and received by about 150 senior executives at the Enron Corporation, over the period 1998 to 2002. Although it is well-known to suffer from various integrity problems, it makes an attractive real data example because it is publicably available and many contemporary readers will be familiar with emailing behaviour. The dataset we analyse was downloaded from http:// bailando.sims.berkeley.edu/enron_email.html and reprocessed for our application. Only emails sent during the year 2001 were retained, because the record appears to be cleanest for that year. Some further effort was then needed to obtain reliable data. For example, many different email addresses can correspond to the same identity since an individual, John Smith say, could appear as any of john.smith, x..smith, jsmith @ either enron.com or ect.enron.com and more.
Following Perry and Wolfe (2013) we discarded emails sent to more than 5 recipients, a subjectively chosen threshold that allows us to focus on inter-personal communications rather than company-wide announcements.
Results will be presented for an individual, hereafter identified as o, who emailed frequently over the year, and for whom there are 12 individuals (of the 150 above) who contact o and that o contacts back. These are referred to by the identifiers 1, . . . , 12.
We shall seek to determine whether i → o triggers o → j, denoted i → o o → j, using only the timing of events. When i = j, a significant test is evidence for reciprocation (or o responding to emails), otherwise it suggests information flow. The point processes generated by i → o and o → j replace A and B respectively in Section 2.

Bayesian null model
To use the results of Section 2, we first require r in (1). In this section we proposes to estimate this function via a generative model, described next. Let D = {t 1 , . . . , t N } denote the times of all of the emails sent by o and D j the times specific to a recipient j. For simplicity D j is modelled as a completely random subset of D, meaning that it can be generated, under the null, by randomly discarding each element of D with a fixed probability. Figure 6 presents o's sent email times, with the y-axis showing the day and x-axis the time of day of each event. This brings out a daily pattern in emailing behaviour; for example between the hours of 1600 and midnight o is markedly less busy, a period of 8 hours where o is presumably away from work (the timezone in which he was operating is not known to us). On the other hand, looking vertically, there some suggestion of behavioural change at a larger time-scale. For example there appears to be a busy period around the middle of year (specifically, July). This visual analysis invites us to separate the likelihood for D as where N k is the number of events in day k and d i = t i mod 24, i = 1, . . . , N are the events collapsed onto a day.
ν 2 was set to 2. An advantage of this approach is a posteriori λ and F are independent, and so posterior samples for each can be obtained in completely separate simulations. In each case these were obtained using reversible-jump Markov Chain Monte Carlo (Green, 1995), using one million samples and retaining every 100th term. (The sample is thinned because the computation below is more expensive than a sampling step.) To obtain an estimate of r, we proceed as follows. We take the averageR of where d(t) is the day on which t occurs, and R n and R d are formed from individual posterior draws of λ and f as  presented. The line isr(t), which we will use in place of r in (1). The model finds a unimodal daily pattern and, for instance, a period of high activity between 20th and the 25th of January. The black circles indicate results retained by FDR set at 10%, the half-circles p-values not retained by FDR but still smaller than 0.05, while the white circles indicate a p-value greater than 0.05.
Next we turn to detecting triggering effects between o's sent and received emails. In Figure 4, the results of testing i → o o → j are shown for i, j = 1, . . . , 12, limiting the range of τ to τ max = 1 week, as discussed in Section 3.1. The black circles are the p-values that are retained for analysis if we use a false discovery rate (FDR) of 10% (Benjamini and Hochberg, 1995). The half-circles are for p-values less than 5% and white are the (non-significant at the 5% level) p-values that remain.
Most of the entries on the bottom-left to top-right diagonal are black, meaning that there is compelling evidence for reciprocation. Because reciprocation is largely to be expected, the two white circles on that diagonal warrant additional inspection: they indicate a lack of evidence for o responding to emails from either 7 or 10.
There is only one email 7 to o and one other from o to 7. They are sent about one month apart (and appear to be unrelated judging by their subject-lines). The p-value is automatically 1 because τ max was set to a week.
10 → o o → 10 This example is more interesting. The p-value is only 0.28 despite there being 14 emails from 10 to o and 9 from o to 10, the most coincidental email times falling in July, about 3 1 2 hours from each other. The reason why no effect is detected is in part because o is estimated to be relatively busy in July, with ρ(July) ≈ 0.13 as opposed to the average 1/12 ≈ 0.08, meaning that we are less sensitive to coincidental timings during that month than at other times. In fact, upon inspecting the subjectlines of 10 → o and o → 10, it does appear as if 10 and o do not reciprocate. For example, the subject-lines of the two most coincidental emails are "FW: Enron Complaint" and "Dunn hearing link?", not obviously related.
Consider now the p-values that were retained by the FDR procedure. Table 1 shows the subject-lines of the emails that the test based its decision on. More precisely, for each of the retained i, j pairs, we find the closest two i → o and o → j events, subject to the former preceding the latter. Thus we have the 'most triggering' email event e, and its subject-line is displayed first. Next, we display the subject-lines of all the emails from o to j that fall within [e, e +τ ], for theτ used by the test. Table 1 shows that the method succeeds in picking out 'real' excitation periods. Consider for example the most significant detection 11 → o o → 11. 4 emails fall withinτ of the 'most triggering' email. These all have the subject-line "RE: DWR -Gas Daily" corresponding quite compellingly to the subject-line of the original email by 11 , "DWR -Gas Daily". Furthermore, the next email in the o → 11 sequence has a different subject-line, "RE: DWR and Edison Meetings", suggesting thatτ fits well.
Although the method largely found evidence of reciprocation, there are some places where real information flow was identified. The email from 9 to o with subject "California Update-Legislative Push Underway" being followed by a string of emails from o to 3 with subject "Re: California Update-Legislative Push Underway" is a particularly compelling example. A concern could be that 3 was simply been 'cc'ed' while o was responding to 9. This is not the case: two of the four emails displayed, the second and the fourth, are sent from o directly to 3 with no other party involved. For reference, there were 28 emails in [e, e +τ ]. Only those with the matching subject-lines are displayed. The vertical dots indicate the position of the omitted emails.
There are admittedly what appear to be false positives in the results, notably the detection 5 → o o → 4. Although there are over 40 hours separating the 'most triggering email' from the two emails that are considered 'triggered', these are the only two emails from o to 4 that occur after the first email from 5 to o, about half a year earlier. The problem is that our model for r does not really capture 'bursty' behaviour -our intuition tells us that the second email, just five hours later, is very likely to be related to the first (and the subject-lines suggest this) rather than adding much evidence to 5 → o triggering o → 4. The step-function that is fit, however, looks for changes on a much larger time scale. Thus we postulate that a more tailored model for r, especially one that captured these local bursts of activity would improve detections.

Conclusion
This article proposes a family of tests for association between two point processes, and analyses their finite-sample and asymptotic properties. Depending on implementation details the procedure can be used to determine different forms of association, including triggering, correlation, inhibition or anti-correlation.
It is remarkable that such a simple and exact technique should exist and not have been discovered until now. Many exciting extensions are now possible. In the analysis of dynamic network data, it may be of interest to find paths or other motifs, rather than simple pairwise dependence. Here the difficulty would be computational as well as statistical, because of the potential complexity of the search. Another important addition would be to incorporate information associated to each point in order to better identify 'pairings' between the events of both processes. This work also leads to an interesting probabilistic question. Consider that the approach relies on the self-superposition and time-transformation of point processes -different tests are obtained depending on the order in which these operations are performed.
To extend results to situations where the null intensity is estimated from the data, the key Table 1: Subject-lines of sent emails that are estimated to be triggered. Further details in main text. obstacle seems to be characterizing point processes that have been time-transformed using an estimated conditional intensity. In particular it would be interesting to investigate conditions under which the resulting points would be more evenly dispersed than Poisson, in some sense to be formalised, and whether this could allow conservative testing.

A On the maximum likelihood estimates
Procedure 1 is established as follows. Firstly, any τ for which K(τ )/n < ρ{T(τ )} can be rejected. To see this, note that for any given τ , (3) is maximised over λ 1 ≥ 0, λ 2 ≥ 0 at λ These estimates must be rejected if λ * 1 < λ * 2 , that is, if K(τ )/n < ρ{T(τ )}. In this case, the maximum over the restriction λ 1 ≥ λ 2 must satisfy λ * 1 = λ * 2 because, while λ * 1 > λ * 2 , we can continue to increase the likelihood by bringing at least one of λ * 1 , λ * 2 closer to K(τ )/n or (n − K(τ ))/n, respectively. But if we can only consider solutions where λ * 1 = λ * 2 , then we are effectively fitting the single parameter model from the null hypothesis. That is automatically worse than choosing τ just large enough that T(τ ) contains all of the events, where we would estimate λ * 1 = n, λ * 2 = 0. The second point is that amongst the τ s that satisfy K(τ )/n ≥ ρ[T(τ )], the likelihood is always improved by shifting τ towards the closest b i − a(b i ) on the left of it and refitting λ 1 , λ 2 , essentially producing a 'sharper fit'. Thus, we can restrict our search on τ to a subset of b 1 − a(b 1 ), . . . , b n − a(b m ) or equivalently restrict our search on ρ{T(τ )} to u 1 , . . . , u n defined in the main text, leading to the solution given.

B Superposition based on conditional intensity
Let L = 1, A = {0, 1/2} and consider a point process B with just one point uniformly distributed on [0, 1). Of course Lemma 1 cannot hold since there can only be one uvariable whereas for λ 1 ≥ λ 2 > 0 the Lemma would admit any number with positive probability. Nevertheless, to illustrate the problem of extending the Lemma to a general point process, consider the distribution of u if r in (1) is replaced by the point process conditional intensity 1/(1 − t), with λ 1 = λ 2 = 1. Let v be a uniform random variable over [0, 1/2). We find that u is distributed as a mixture − log(1 − v) (an exponential variable conditioned on being smaller than log 2) with probability 1/2 and − log(1 − v) 1−2v (a non-standard distribution).