Large-scale multivariate point process data are fast emerging in a wide range of scientific and business applications. Learning the latent network structure from such data has become an increasingly important task. For instance, one may wish to estimate the neuronal functional connectivity network based on spiking times (i.e., times when a neuron fires) recorded from ac o l l e c t i o no fn e u r o n s( F a r a j t a b a re ta l .2015), or to estimate the financial network based on trading times recorded for a collection of stocks (Linderman and Adams 2014). Both the neuron spiking times and the trading times can be viewed as realizations from multivariate point processes. To characterize the latent interactions between the different point processes, a useful class of models is the multivariate Hawkes process (Hawkes 1971). The multivariate Hawkes process is a mutuallyexciting point process, in which the arrival of one event in one point process may trigger those of future events across the different processes. Because of its flexibility and interpretability, the multivariate Hawkes process has been widely used in many applications, such as social studies (Zhou, Zha, and Song 2013), criminology (Linderman and Adams 2014), finance (Bacry et al. 2013) and neuroscience (Okatan, Wilson, and Brown 2005). In the network setting, each component point process of the multivariate Hawkes process is viewed as a node, with a directed edge connecting two nodes indicating an event in the source point process increases the p r o b a b i l i t yo fo c c u r r e n c eo ff u t u r ee v e n t si nt h et a r g e tp o i n t process.
Despite the popularity of the multivariate Hawkes process, there is a need of new statistical theory and methodology for its broader applications. First, most existing theoretical results for the Hawkes process are derived using a cluster process representation of the process. This cluster process representation by its definition depends on the mutually excitation a s s u m p t i o n ,t h a ti s ,t h ea r r i v a lo fo n ee v e n ti n c r e a s e st h e p r obab ili tyo foccurr enceo ffu t ur eeven ts(H a wk esa ndOak es 1974; Hansen, Reynaud-Bouret, and Rivoirard 2015). Howe v e r ,s u c ha na s s u m p t i o nm a yn o tb ev a l i di nc e r t a i na p p l ications. For example, it is well known that the firing activity of one neuron can inhibit the activities of other neurons (Amari 1977). A more flexible model should allow both excitatory and inhibitory effects, which renders the cluster process representation infeasible. Second, most existing models assume that the background intensities, that is, the baseline arrival rates of events from the different component processes, are constant in time. Under this assumption, the multivariate Hawkes process satisfies a stationary condition (e.g., Brémaud and Massoulié 1996). However, assuming constant background intensities may also be too restrictive. For example, stock trading activities tend to be much higher during market opens and closes (Engle and Russell 1998), and the associated background intensities are therefore not constant in time. A multivariate Hawkes process with constant background intensities may not fit such data well (Chen and Hall 2013). A more flexible approach instead should allow the background intensities to be time-varying. For such nonstationary models, new development on both theory and methodology is needed, as most existing results are established assuming the underlying process to be stationary.
Some existing work have considered broadening the class of Hawkes process models. Specifically, Brémaud and Massoulié (1996), Chen et al. (2017), and Costa et al. (2018) considered a class of nonlinear Hawkes processes that allows both excitatory and inhibitory effects. A thinning process representation was used to investigate the properties of the proposed process. However, these work focused on processes with constant background intensities and the thinning representation technique depended critically on the stationarity condition. Some recent work also considered nonstationary Hawkes processes. Lewis and Mohler (2011), Chen and Hall (2013), and Roueff, Von Sachs, and Sansonnet (2016) considered Hawkes processes with time varying background intensities. However, they only considered univariate processes, and only with excitatory effects. Lemonnier and Vayat is (2014) considered a multivariate Hawkes process with time varying background intensities. However, they focused on an approximate optimization algorithm for model estimation, and did not provide any theoretical results.
In this article, we propose a flexible class of multivariate Hawkes process that admits time-varying background intensitiesandallowsbothexcitatoryandinhibitoryeffects.W eshow the existence of a thinning process representation of this nonstationary process. Such a result has not yet been established in the literature, and it enables our subsequent theoretical analysis. To estimate the network structure, we consider a computationally efficient penalized least squares estimation, in which both the background intensities and the transfer functions are approximated using basis functions. We establish theoretical properties of the penalized least squares estimator in the high-dimensional regime, where the dimension of the multivariate process p can grown much faster than the length of the observation window T. Specifically, we investigate the following properties in our analysis:
1. (Concentration inequalities.) We establish concentration inequalities for the first and second order statistics of the proposed Hawkes process. Such inequalities are established using a new thinning process representation result for nonlinear and nonstationary Hawkes processes. 2. (Non asymptotic error bound.) Under certain regularity conditions, we establish, in the high-dimensional regime, the non asymptotic error bound of the intensity functions estimated using basis approximations. Specifically, we verify that the design matrix satisfy a restricted eigenvalue condition and a bounded eigenvalue condition for the diagonal blocks; these bounds on eigenvalues depend on the number of basis functions. 3. (Network recovery.) We show that, under certain regularity conditions, our proposed estimation method can consistently identify the true edges in the network with probability tending to one. Moreover, we propose a consistent generalized information criterion (GIC) for regularizing parameter selection. 4.(T estforbackgroundintensity .)W eproposealeastsquares based statistic for testing if the background intensity is constant in time. Specifically, we show that the null distribution of the test statistic is asymptotically χ 2 and the test is powerful against alternatives.
It is worth mentioning that there is another class of approaches that estimate the latent network structure from high dimensional multivariate point process data (Zhang et al. 2016;Vinci et al. 2016Vinci et al. , 2018)). These methods divide the observation windowintoanumberofbins,andmodelthenumberofevents in each bin. The network structure is estimated using methods such as correlation of event counts (Vinci et al. 2016), regularized generalized linear models (Zhang et al. 2016), or Gaussian graphical models (Vinci et al. 2018). The heuristic binning procedure may result in information loss, for example, short-term excitatory effects may be overlooked if the bins are chosen to be too wide. Choosing an appropriate binning procedure remains a challenging task.
The rest of the article is organized as follows. Section 2 i n t r o d u c e st h ep r o p o s e dm o d e l ,a n dSection 3 describes the model estimation and selection. The aforementioned theoretical results are detailed in Sections 4. Section 5 includes simulation studies. The detailed analysis of a neuron spike train dataset is presented in Section 6. A short discussion section concludes the article.
when there is no ambiguity). Let f (k) denote the kth derivative of a function f when such a derivative exists. For a matrix A ∈ R m×n ,w eu s e A 2 , A max and A ∞ to denote its spectral norm, maximum entry-wise ℓ 1 norm and maximum row-wise ℓ 1 norm, respectively. We write
[n]={ 1, 2, ..., n} and let x denote the largest integer less than x.F o ras e tS, we use |S| to denote its cardinality. We write 1 n to denote a length-n vector of 1, I n×n to denote a n × n identity matrix, diag{d 1 , ..., d n } to denote a n×n diagonal matrix with diagonal elements d 1 , ..., d n ,anduseσ min (•) and σ max (•) to denote the smallest and largest eigenvalues of a matrix, respectively. For two positive sequences a n and b n ,w ri t ea n
Consider a directed network with p nodes. For each node j ∈[p], we observe its event locations {t j,1 , t j,2 , ...} in the time interval [0, T] such that 0 < t j,1 < t j,2
< ••• ≤ T.F o rn o d ej,l e t the associated counting process be N j (t) =| { i : t j,i ≤ t}|, t ∈[ 0, T].W r i t eN = (N j ) j∈[p] as the p-variate counting process. Let H t denote the entire history of N up to time t,and write N j ([t, t + dt)) as dN j (t).Thep-variate intensity function
We propose a flexible class of Hawkes processes with intensity functions defined as
where h(•) : R → R + is a link function and it is assumed to be θ -Lipschitz (see Assumption 1), ν j (•) : R + → R + is the timevarying background (or baseline) intensity function of the jth process, and ω j,k (•) : R + → R is the transfer function that characterizes the effect of the kth process on the jth process. Specifically, i. ω j,k (s)>0 corresponds to excitatory effect, that is, an event in process k increases the probability of event occurrence in process j at a time distance of s. ii. ω j,k (s)<0 corresponds to inhibitory effect, that is, an event in process k decreases the probability of event occurrence in process j at a time distance of s. iii. ω j,k (s) = 0 corresponds to no effect, that is, an event in process k has no effect on the event occurrence in process j at a time distance of s.
The proposed model in (1) considers a time dependent background intensity function instead of the constant background intensity considered in existing multivariate Hawkes process models (Chen and Hall 2013;Hansen, Reynaud-Bouret, and Rivoirard 2015;Bacry et al. 2020;Wang, Kolar, and Shojaie 2020). Consequently, the proposed Hawkes process model is nonstationary, and its analysis requires new theoretical tools, which will be introduced in Section 4.
Let the directed network G(V, E) summarize the relationships between the p component processes. Specifically, let V = {1, 2, ..., p} be the set of p nodes and E be the set of edges such that
where ω j,k 's are the transfer functions in (1). Therefore, (j, k) ∈ E if and only if the kth process has an excitatory or inhibitory effect on the jth process.
Next, we introduce a set of regularity conditions on the background intensities and transfer functions in (1).
Assumption 1. Let be a p × p matrix with jk = ∞ 0 |ω j,k (t)|dt, j, k ∈[ p] and assume that σ max ( ) ≤ σ < 1. Moreover, assume that h(•) is a θ -Lipschitz link function with θ ≤ 1, and the background intensity functions are bounded, that is, 0 <ν j (t) ≤ ν, j ∈[p], for some positive constant ν.
Assuming the Lipschitz constant θ to satisfy θ ≤ 1i sn o t restrictive. For example, if h(•) is K 0 -Lipschitz for some K 0 > 1, we can reparameterize (1)b ysetting h(x) = h(x/K 0 ), ν j (t) = K 0 ν j (t) and ω j,k (t) = K 0 ω j,k (t). In this reparameterized model, h(•) is 1-Lipschitz. The Lipschitz condition on the link function was also considered in Massoulié (1998) and Chen et al. (2017). Assumption 1 implies that h{ν j (t)} is bounded as Lipschitz functionsareboundedonboundedsupports.
Define the mean intensity of (1)a s λj (t) = E{dN j (t)}/dt. Under Assumption 1,themeanintensity λj (t) is upper bounded. This can be shown in three steps. First, define a p-dimensional Hawkes process N * = (N * j ) j∈ [p] with intensity function
where ν * is a positive constant such that ν * ≥ h{ν j (t)} for any j and t,andω j,k 'sareasdefinedin(1). By Brémaud and Massoulié (1996), the point process defined in (2) satisfies a stationary condition under Assumption 1. Next, write the mean intensity of (2)a s * = ( * 1 , ..., * p ) ,w h e r e * j = E{dN * j (t)}/dt. Correspondingly, we have
where ν * = (ν * , ..., ν * ) ∈ R p and ω(t) ∈ R p×p ,w i t h {ω(t)} jk = ω j,k (t).Themeanintensity * in (3)canberewritten as * = ∞ k=0 k ν * ,whichisupperboundedgivenσ < 1in Assumption 1.Finally ,itcanbeshownthatthemeanintensity of (1), that is, λj (t), is upper bounded by * j (see Lemma S6 and its proof in the supplementary materials). Consequently, λj (t) is also upper bounded under Assumption 1.
From the observed event locations in [0, T],ourobjectiveisto estimate the intensity functions λ j (t), j ∈[ p].Furthermore,by identifying the nonzero transfer functions ω j,k 's in the estimated intensity functions, we can estimate the structure of the directed network G(V, E). In this section, we consider the link function to be h(x) = max(0, x) in (1). It is seen that this link function is a 1-Lipschitz function. To ease notation, we write
To estimate the intensity functions, one may consider a likelihood function based approach (Ogata 1981;C h e na n dH a l l 2013; Zhou, Zha, and Song 2013). However, minimizing the negative log-likelihood function may require a very involved and computationally intensive iterative procedure (Veen and Schoenberg 2008). To improve the estimation efficiency, we consider a least squares loss based estimation approach. That is, we consider the following loss function
The least squares loss comes from the empirical risk minimization principle (van de Geer 2000)a n dh a sb e e nf a i r l yc o mmonly considered in estimating point process models (Hansen, Reynaud-Bouret, and Rivoirard 2015;C h e ne ta l .2017;B a c ry et al. 2020). We later show that (5)c a nb es e p a r a t e di n t op objective functions that can be estimated individually, which significantly reduces the computation cost.
We consider a nonparametric estimation of ν j (t) and ω j,k (t), j, k ∈[ p], using B-spline approximations. Given [a 1 , a 2 ]⊂R and a set of
),letS K,l be the space of polynomial splines of degree l ≥ 1 consisting of functions satisfying: (i) restricting to each interval
the function is l times continuously differentiable (Stone 1985). Such a space S K,l is of dimension m = K + l (Schumaker 2007)a n da ss u c h ,l e t{φ 1 (t), ..., φ m (t)} be the normalized B-spline basis of S K,l .W h e nl = 1, the basis is a set of K + 1 step functions with jumps at knots (Stone 1985). In our procedure, we approximate the background intensity ν j (t) with an m 0 -dimensional normalized B-spline basis φ 0 (t) = (φ 0,1 (t), ..., φ 0,m 0 (t)) ,s u c ht h a tν j (t) = β j,0 φ 0 (t) + r j,0 (t), where β j,0 ∈ R m 0 and r j,0 (•) denotes the approximation residual. Furthermore, we approximate the transfer functions ω j,k (t)
where β j,k ∈ R m 1 and r j,k (•) denotes the approximation residual. The dimensions and degrees of the bases φ 0 (t) and φ 1 (t) are allowed to be different for more flexibility in characterizing the background intensities and transfer functions. For example, one may use cubic B-splines to approximate the background intensities and step functions to approximate the transfer functions (Hansen, Reynaud-Bouret, and Rivoirard 2015). The choices for the number and locations of knots are discussed in Section 3.1.
Write β j = (β j,0 , β j,1 , ..., β j,p ) . We define α j = (α (j,0) , α (j,1) , ..., α (j,p) ) such that α (j,0) ∈ R m 0 with
and
where the component G (k 1 ,k 2 ) is defined as
With the above expressions for β j , α j and G, we define
Some straightforward algebra shows that the loss function in (5)c a nb ew r i t t e na s p j=1 ℓ j (β j ). We note that both α j and G arecalculatedbasedontheobservedeventlocationsandthe prespecified basis functions, that is, φ 0 (t) and φ 1 (t). Therefore, to estimate the background intensities and transfer functions, we can directly optimize (7)w i t hr e s p e c tt oβ j .S i n c et h el o s s function
w h e r e ω j,k 's are t he estimated transfer functions. Note that if ω j,k = 0, then all coefficients associated with ω j,k are zero (i.e., β j,k = 0). Thus, to encourage sparsity in the estimated network, we impose a standardized group lasso penalty on β j with the coefficients in β j,k grouped together, k ∈[ p]. Specifically, we consider the following optimization problem min
The penalty term
is an extension of the standardized group lasso penalty (Simon et al. 2013). This optimization problem in (8)i sc o n v e xa n dc a nb ee ffi c i e n t l y solved using a block coordinate descent algorithm (Simon et al. 2013). The terms α 1 , ..., α p ,a n dG c a nb ec o m p u t e du s i n g standard numerical integration methods and such calculations c a nb ec a r r i e do u tb e f o r ei m p l e m e n t i n gt h eb l o c kc o o r d i n a t e descent algorithm.
Our proposed estimation procedure involves a number of tuning parameters, including the numbers of B-splines (i.e., m 0 and m 1 ) for approximating the background and transfer functions, respectively, knots locations for the B-splines, and tuning parameter η j 's in the penalized least squares estimation in (8). Cross-validation procedures may not be appropriate for tuning parameter selections under our setting as the proposed process is nonstationary due to the time-varying background intensity in (1). As such, the data cannot be divided into training and validation sets in a straightforward manner.
Given m 0 and m 1 , we let the knots be evenly distributed (Ravikumar et al. 2009;H ua n g,H o r o wi t z ,a n dW e i2010). For m 0 and m 1 , theoretical conditions in Theorem 3 can guide their empirical choices. In Section B1, we describe a heuristic procedure for selecting m 0 and m 1 ; a similar procedure was considered in Kozbur (2020). In Section B2, we show that this heuristic procedure achieves good performance; additionally, we demonstrate that the estimation accuracy is not overly sensitive to the choices of m 0 and m 1 .O n c eφ 0 (t) and φ 1 (t) are determined, we then move to select η j 's.
The tuning parameter η j in (8) controls the sparsity of β j , which in turn controls the sparsity of the estimated network. To select η j , we propose a generalized information criterion (GIC) defined as
where ℓ j (•) i sa sdefin edin(7), κ j = T/N j {(0, T]} is a scaling parameter, β j is estimated from (8)w i t hη j ,a n dα T > 0i s a parameter that scales with T and p.A sℓ j ( β j ) is the squares loss and not the log-likelihood function, the GIC is not directly comparable to the likelihood based selection criteria such as the BIC or extended BIC (Schwarz 1978;C h e na n dC h e n2008).
In Theorem 5, we show that the proposed GIC is consistent given appropriate choices of α T ,s u c ha sO((log p) 2 log T).I n Section 5,weevaluatetheefficacyoftheproposedGICandshow it achieves satisfactory performance.
In this section, we first show the existence of a thinning process representation of the proposed nonlinear and nonstationary Hawkes process. We then establish concentration inequalities f o rt h efi r s ta n ds e c o n do r d e rs t a t i s t i c so ft h ep r o p o s e dp o i n t process. These results are useful in the subsequent analysis of the estimated intensity functions. Next, we establish the nonasymptotic error bound of the estimated intensity functions and show that our method can consistently identify the true edges in the network. Lastly, we propose a test statistic for testing if the background intensities are constant in time. We derive its asymptotic null distribution and show the test is powerful against alternatives. All proofs are collected in the supplementary materials.
Many existing theoretical analyses rely on the cluster process representation of the Hawkes process (Hawkes and Oakes 1974;Hansen, Reynaud-Bouret, and Rivoirard 2015;Bacry et al. 2020), which needs the transfer functions to be nonnegative. Brémaud and Massoulié (1996) employed a thinning process representation of the Hawkes process that permitted negative transfer functions, but required a stationarity condition. As such, the results in Brémaud and Massoulié (1996)a r en o t directly applicable to our problem. Next, we first show the existence of a thinning process representation of the proposed nonlinear and nonstationary Hawkes process in (1).
Let N = (N j ) j∈[p] be a p-variate homogeneous Poisson process on R 2 with intensity 1. Let λ (0)
and
as follows:
where h(•), ν j and w j,k are as defined in (1), and N j ([0, λ (n+1) j (t)]×dt) denotes the number of points for
is the intensity function of the point process N (n) (t). Next, we show that the sequence {N (n) } ∞ n=1 in (10) converges in distribution to the Hawkes process N with intensity function (1).
Theorem 1. Let λ(t) be as defined in (1)satisfyingAssumption 1. Let {λ (n) (t)} ∞ n=1 and {N (n) } ∞ n=1 be sequences as defined in (10). Then, it holds that (a) λ (n) (t) converges to λ(t) almost surely for any t, (b) N (n) ∞ n=1 converges in distribution to N with intensity (1).
Theorem 1 shows the existence of a thinning process representation of the proposed nonstationary Hawkes Process. It provides a theoretical guarantee, analogous to Massoulié (1998), for nonstationary multivariate Hawkes processes. This new result is critical in our subsequent theoretical analysis.
Assumption 2. Assume that there exists λ max > 0s u c ht h a t λ j (t) ≤ λ max for any t and j and ω j,k , j, k ∈[ p] are bounded functions with a bounded support [0, b] for some b > 0.
This condition first assumes that the intensities are upper bounded. One example of such processes is when the link function h(•) is upper bounded by a positive constant; see also Section 7. Assumption 2 also assumes that the transfer functions ω j,k 's have a bounded support. The bounded support assumption has been fairly commonly considered in the analysis of multivariate Hawkes process (Hansen, Reynaud-Bouret, and Rivoirard 2015;Costaetal.2018
Assumption 3. There exists ρ ∈ (0, 1) such that
This assumption requires that has bounded column sums, which prevents the intensity function from concentrating on any single process.
Recall that H t denotes the history of N up to time t.ForH tpredictable functions f 1 (•) and f 2 (•), define
Theorem 2. Consider a Hawkes process on [0, T] with intensity as defined in (1)s a t i s f y i n gAssumptions 1-3.L e tf 1 (t) be a bounded function and f 2 (t) be a bounded function on a bounded support. Then, for k ∈[p],itholdsthat
where c 1 , c 2 ,andc 3 are positive constants. For any j, k ∈[ p],it holds that
where c 1 , c 2 ,andc 3 are positive constants.
The proof of Theorem 2 is provided in the supplementary materials. Our proof strategy for the concentration results in Theorem 2 follows from that in Chen et al. (2017). Specifically, as in Chen et al. (2017), we first define a coupling process of N, which is used to bound the temporal dependence of N.Th en, a Bernstein type inequality for weakly dependent sequences (Merlevède, Peligrad, and Rio 2011)isusedtoobtainthedesired results. The main difference in the proof is that the validity of the thinning process representation in Chen et al. (2017)isensured by Massoulié (1998), established under a stationarity condition. In the nonstationary case, the validity of the thinning process representation is established in Theorem 1.
We remark that the concentration results for stationary processes in Chen et al. ( 2017)m a yn o tb ed i r e c t l ya p p l i c a b l et o establish ( 11)-( 12) under the nonstationary case we considered, even though N is dominated by a stationary process (see Lemma S6), as no existing theoretical results, to our knowledge, establish the concentration inequalities of a target process using directly the concentration inequalities of its dominating process. Theorem 2 is useful in the ensuing theoretical analysis that derives the non-asymptotic error bound of the estimated intensity functions and establishes edge selection consistency. The next corollary is a direct consequence of Theorem 2.
Corollary 1. Consider a Hawkes process on [0, T] with intensity a sd e fi n edi n( 1)sa ti s fy i n gAssumptions 1-3. Considering the matrix G defined in (6), we have
where c 4 , c 5 ,andc 6 are positive constants.
The result in Corollary 1 is a direct consequence of Theorem 2, once we show that the entries in G are first and second order statistics of the proposed Hawkes process.
In this section, we derive the non asymptotic error bound of the estimated intensity function in the diverging p regime. @To simplify notation, we define (t) = ( 0 (t), 1 (t), ..., p (t)) , where 0 (t) = φ 0 (t) and k (t)
Recall the first order mean intensity function λk (u) is defined as λk
Denote the p × p covariance function as C 0 (u 1 , u 2 ),suchthat,
entry is defined as
When Hawkes 1971). Thus, the complete covariance matrixcanbewrittenas (Hawkes 1971). Note that C(u 1 , u 2 ) is in general not symmetric (Li and Zhang 2011). Specifically, C(u 1 , u 2 ) is symmetric when u 1 = u 2 ;w h e nu 1 = u 2 , C(u 1 , u 2 ) i st h ec r o s s -c o v a r i a n c ef u n c t i o n ,w h i c hi sn o t necessarily symmetric, and it holds by definition that
Assumption 4. Assume there exist constants min , max > 0 such that, λk (t) ≥ min and λ(2)
This condition assumes that the first and second order mean intensities are bounded. The lower-bounded condition on the mean intensity λj (t) c a nb es a t i s fi e dw h e nt h ei n h i b i t o r y effect from the negative transfer functions is not excessive when compared to the background intensity and the excitatory effect from the positive transfer functions. The nonnegative definite assumption of C 0 (u 1 , u 2 ) holds true for many commonly used univariate point process models (Guan, Jalilian, and Waagepetersen 2013). In the stationary multivariate Hawkes process case, Bacry and Muzy (2016)showedthatC 0 (u 1 , u 2 ) is directly related to the solution to an integral equation involving the transfer functions; the integral equation can be numerically solved and an estimate of C 0 (u 1 , u 2 ) can therefore be obtained. In our nonstationary multivariate Hawkes process setup, C 0 (u 1 , u 2 ) may instead be estimated through parametric bootstrap (see Section B4). Validity of the nonnegative definite assumption of C 0 (u 1 , u 2 ) canbesubsequentlyassessedusingan estimated C 0 (u 1 , u 2 ).
Assumption 5. Assume that there exist β j = ( β j,0 , β j,1 , ..., β j,p ) ∈ R m 0 +pm 1 , j ∈[ p], and a smoothness parameter d ≥ 2suchthat,forsomepositiveconstantsC 1 , C 2 and C 3 ,
with probability at least 1 -C 2 pT exp(-C 3 T 1/5 ),wherem 1 m 0 and β j,k = 0 for k / ∈ E j .
This condition assumes that the true intensity function can be well approximated by the basis functions, in that residuals from the truncated basis approximation decreases at a polynomial rate of the number of basis functions. The d ≥ 2i s a smoothness parameter for the background intensities ν j (t)'s and transfer functions ω j,k (t)'s. While this parameter may differ between ν j (t)'s and ω j,k (t)' s ,i ti sa s s u m e dt ob et h es a m et o simplify notations in our analysis. Condition (15)canbeverified when, for example, h(x) = x and the approximation errors
,wh e r eb is as defined in Assumption 2; see a detailed proof of this statement in Section A10 and also Section 7.Such approximation errors hold for B-spline basis (Stone 1985)o r trigonometric basis (Tsybakov 2008) when the target functions belong to certain function classes. For example, when ω j,k is dsmooth (Chen 2007), that is, |ω (l) j,k (t)ω (l) j,k (s)|≤c|t -s| d-l , where l = d and c is some positive constant, there exists b 1985). We refer to C hen (2007)a n dT s y b a k o v( 2008) for thorough reviews of basis approximations and truncation errors.
Next we establish the non-asymptotic error bound of the estimated intensity functions.
Theorem 3. Consider a Hawkes process on [0, T] with intensity a sd e fi n e di n( 1)s a t i s f y i n gAssumptions 1-5.F o rj ∈[ p], let λ j (t) = (t) β j ,w h e r e β j is estimated from (8). Given η j = (C 2 log p/T) 1/2 , s = o(T 2/5 ),l o gp = O(T 1/5 ) and sm 1 = O(T 4/5 ),wehave,forj ∈[p],
holds with probability at least 1-C 3 p -2 -C 4 p 2 T exp(-C 5 T 1/5 ), where C 2 , C 3 , C 4 ,a n dC 5 are positive constants, and C 1 is as definedin(15).
T h ee r r o rb o u n do nt h er i g h th a n ds i d eo f( 16)c o n s i s t so f two terms. The first term comes from the B-spline basis approximation error (i.e., bias from approximating the nonparametric background and transfer functions using basis functions) and the second term comes from the statistical error (i.e., stochastic error in estimating the intensity functions). It is seen that when sT/ log p = o(m 2d 1 ),t h eb i a st e r mC 1 (s + 1) 2 m -2d 1 would become negligible when compared to the statistical error term. When, for example, s = O(1), d = 2andm 1 T 1/5 ,theerror bound in ( 16)r e d u c e st oO(T -4/5 + log p/T),wh i c hi sc o mparable with the estimation error in sparse additive regressions (Raskutti, Wainwright, and Yu 2012).
Two key ingredients in the proof of Theorem 3 are establishing an upper and lower bounded eigenvalue condition for G (k,k) (see Lemma S10) employed in the standardized group lasso penalty in (8) and a restricted eigenvalue condition for G under the group lasso setting (see Lemma S11). Establishing these two conditions under the proposed nonstationary process is nontrivial; it requires a delicate analysis that combines properties of the basis functions and concentration inequalities of first and second order statistics of the proposed process. Combining these two ingredients and a martingale central limit theorem for counting processes (van de Geer 1995), we are able to derive the result in Theorem 3.W en o t et h a ti ft h eb a s i sa p p r o x i m a t i o n error condition in Assumption 5 is not satisfied, we may replace the first term in the error bound ( 16), that is,
and Theorem 3 holds with an error bound of 32(R m 0 ,m 1 + 9sλ max log p/T).
In this section, we show that, under certain regularity conditions, our proposed method can consistently identify the true edges in the network with probability tending to one. Recalling
where
a n d γ min and γ max are constants as defined in Lemma S10. This is the irrepresentable condition (Zhao and Yu 2006) under our setting and it is a condition on covariances between the component processes. Considering the jth component process, this condition stipulates that the k (t) from non neighbors of j (i.e., k / ∈ E j )h a ss m a l lc o v a r i a n c e sw i t h k (t) from neighbors of j (i.e., k ∈
and some constant c 0 > 0; see a detailed proof in Section A11. The condition canbefurtherrelaxediftheada ptivelasso(Zou2006;H uang, Horowitz, and Wei 2010) penalty term is considered and we plan to investigate this extension in our future work. The next condition is a minimal signal condition.
Assumption 7. There exists a constant β min > 0s u c ht h a t β j,k 2 ≥ β min for k ∈ E j ,where β j is as defined in Assumption 5.
Note that this condition is not placed on β j,0 since β j,0 is not included in the penalty term. This minimal signal condition can be relaxed; see discussions in Section 7.
Theorem 4. Consider a Hawkes process on [0, T] with intensity asdefinedin(1)satisfyingAssumptions 1-7. Assume that η j =
with probability at least 1 -2C 3 p -2 -3C 4 p 2 T exp(-C 5 T 1/5 ), where C 2 , C 3 , C 4 ,andC 5 are as in Theorem 3. This result establishes selection consistency. The condition s 2 T/ log p = O(m 2d 1 ) is needed in selection to ensure the bias term O(s 2 m -2d 1 ) does not dominate the group-wise estimation error O(log p/T). The conditions in Theorem 4 are satisfied when, for example, s = O(1), d = 2, log p T 1/5 and m 1
T 1/5 .W eno tetha tselectio nco n sistencywasalsost udi e di nC h e ne ta l .( 2017), under a stationary Hawkes process setting. In comparison, our result is established under the more flexible nonstationary setting. Moreover, our result significantly relaxes a restrictive condition in Chen et al. (2017). Specifically, Assumption 7 (second equation) in Chen et al. (2017), after some simplification, would require T to be upper bounded. This result on the selection consistency has an important implication in practice, as it ensures that our method can correctly identify the true edges in the latent network.
Next,weinvestigatetheselectionconsistencyoftheproposed GIC in (9). We use E η j j to denote the estimated E j with tuning parameter η j .Letη max and η min be, respectively, the upper and lower limits of the tuning parameter η j ,w h e r eη max can be easily chosen such that E η max j is empty and η min can be chosen such that E η min j is sparse, and the corresponding model size
| satisfies conditions in Theorem 5. We partition the interval [η min , η max ] into two subsets
corresponding to η j 's that result in under-fitted and over-fitted models, respectively. The next result states that the proposed GIC is consistent in model selection.
Theorem 5. Consider a Hawkes process on [0, T] with intensity a sd e fi n e di n( 1)s a t i s f y i n gAssumptions 1-5 and 7 and that
The assumption that there exists η * j ∈[η min , η max ] such that E η * j j = E j is satisfied, for example, by the result in Theorem 4, which would additionally require Assumption 6.T h i si st r u e by noting | E η max j |=0, | E η min j |=s 0 , s = o(s 0 ) and the size of the selected model decreases as η j increases (Zhang, Li, and Ts ai 2010). The main challenge in establishing Theorem 5 is the large number of candidate models in the over-fitted case, which increases combinatorially fast with p.T oo v e r c o m et h i s challenge, we introduce a proxy criterion on a support of size s 0 < p (Zhang, Li, and Tsai 2010); see proof details in Section A7. The two conditions on α T specify a range that ensures consistency. Specifically, s log p/α T = o(1) suggests that α T should diverge adequately fast such that the true model is not dominated by over-fitted models. On the other hand, sm 1 α T /T = o(1) restricts the rate of divergence for α T based on the size of thetruemodels and observation window length T.Ifwetake, for example, α T (log p) 2 log T, both conditions on α T are met. While other choices of α T can also satisfy both conditions, we have chosen a uniform choice α T (log p) 2 log T in our empirical investigations. Moreover, we take s 0 T 2/5 / log log T and choose η max such that that | E
In this section, we consider the problem of testing if the background intensities of the proposed Hawkes process are constant in time. If the background intensities ν j (t)'s in (1)areconstant, under Assumption 1, a stationary process whose intensity follows (1) exists (Brémaud and Massoulié 1996).
In our model, the background intensity ν j (t) for the jth processisrepresentedasφ 0 (t)β j,0 .Withoutlossofgenerality,let the first term in the basis φ 0 (t),thatis,φ 01 (t),betheconstant term. Testing if ν j (t) isconstantintimecanthenbeformulated as testing the following hypotheses:
where
∈ R m 0 ×m 0 .T h et e s t in (17) can detect any fixed departure in ν j (t) from a constant provided that m 0 is sufficiently large (Fan, Zhang, and Zhang 2001). Recall that β j is obtained from
where ℓ j (β j ) is defined as in (7). Next, letting supp 1 (β j ) ={k ∈
[p] : β j,k = 0}, and define the refitted estimator β 1 j and the restricted estimator β H 0 j under H 0 as
Note that the number of basis functions used in ( 18)and( 19) may not be the same; see discussion after Theorem 6.F i n a l l y , we define the test statistic as
where ℓ j (β j ) is defined as in ( 7). The following theorem states the asymptotic null distribution of the test statistic.
Theorem 6. Assume that all conditions in Theorem 4 are satisfied and the m 1 used in ( 19)satisfys 3 T = o(m 2d+1 1
). Under H 0 , we have that
where λj = E{dN j (t)}/dt is a constant under H 0 .
Theorem 6 is derived for a two-step procedure, where the first step involves a regularized estimation in (18)andrequires appropriate regularization to achieve selection consistency and t h es e c o n ds t e pi n v o l v e sr e fi t t i n gb a s e do nt h es e l e c t e ds e t of edges from step 1. The selection consistency in step 1 is ensured by Theorems 4-5, and the tuning parameter η j used to calculate β j is selected following the proposed GIC in (9). In step 2, the condition
) is needed such that the bias from approximating the background and transfer functions using basis functions is asymptotically negligible relative to variance of the test statistic. Such a condition is usually referred to as under-smoothing (Chen 2007)a n di sc o m m o n in nonparametric regression testing problems. For instance, if m 1 T 1/5 is used in (18)(seedi scu s sio no fTheorem 4), the under-smoothing condition s 3 T = o(m 2d+1 1 ) will be satisfied if we multiply m 1 by, for example, a factor of T 1/20 .W edonot assume m 0 is fixed when estimating β 1 j ,thatis,m 0 may increase with T.Asm 0 increases, we may alternatively write the limiting distribution as (S j / λj -m 0 +1)/(2m 0 -2) 1/2 → d N (0, 1),which does not depend on m 0 . In practice, given the m 0 and m 1 used in network structure estimation (see Section B1), we multiply both of them by an under-smoothing factor (e.g., T 1/20 ), which results in the under-smoothed m 0 and m 1 used in the testing procedure.
Based on the result in Theorem 6,wewo uldr ejectthen ull H 0 : Aβ j,0 = 0 if S j / λj ≥ z 1-α ,wher ez α is the αth quantile of χ 2 m 0 -1 and λj is estimated with λj = 1 T T 0 λ j (t)dt.A st h e limiting distribution in Theorem 6 is derived under perfect model recovery, ensured by Assumptions 6-7,cautionsshould be exercised when performing this test, as the testing procedure may be invalid when perfect model recovery does not hold. We plan to relax this assumption on perfect model recovery in our future research; see discussions in Section 7.
Next, we discuss the asymptotic power of our proposed test. The following theorem provides a lower bound on the growth rate of the test statistic under alternatives.
Proposition 1. Assume that all conditions in Theorem 4 are satisfied. For any alternative H 1 such that A β j,0 2 /(s 2 m 1 log p/T) 1/2 →∞,wehave
This result shows that the growth rate of S j under the alternative is at least s log p, while p diverges. The asymptotic null distribution in Theorem 6 and the growth rate under the alternative together suggest that the null and the alternative hypotheses are well separated, and our proposed test is asymptotically powerful against alternatives. Moreover, the test is locally powerful, that is, A β j,0 2 isallowedtotendto0aslongasitdecreasesnofaster than the rate of (s 2 m 1 log p/T) 1/2 .
Compared with the inferential result in Chen and Hall (2013) f o rt h em a x i m u ml i k e l i h o o de s t i m a t o ro fu n i v a r i a t eH a w k e s process, Theorem 6 is derived for least squares based estimation. Hence, the specific form of the martingale, the characterization of its large jumps and the asymptotic variance are different when applying the martingale central limit theorem. Moreover, as S j calculates the difference between two least squares losses, a careful treatment was needed to derive its asymptotic expansion; see Step 1 in the proof. Finally, Chen and Hall (2013)c o nsider parametric forms of the background intensity and transfer function while we consider a nonparametric estimation using B-splines. Consequently, our analysis is challenged by a bias term of T 0 Ēj (t) λ j (t)λ j (t) dt in the test statistic, which is bounded by combining properties of the B-spline basis and resultsinLemmasS2andS11.
In this section, we carry out simulation studies to investigate the finite sample performance of our proposed method, and to compare with existing solutions. We consider three simulation settings. In Simulation 1, we investigate the estimation accuracy of our proposed method; in Simulation 2, we investigate the network edge selection accuracy; in Simulation 3, we evaluate t h es i z ea n dp o w e ro fo u rp r o p o s e dt e s to fh y p o t h e s i s .D u e to space limitations, results from Simulation 3 is relegated to Section B3 of the supplement. We refer to our proposed method for nonstationary Hawkes processes as NStaHawkes.InSimulations 1 and 2, we compare our method with Chen et al. all transfer functions are positive, corresponding to excitatory effects. The second setting considers the network in Figure 1(b), with both positive and negative transfer functions, corresponding to excitatory and inhibitory effects, respectively. Let the network edge set be E ={ (k,1), k = 2, ...,11}.Thebackground intensity functions and transfer functions for each setting are as follows:
Setting 1.1: ν 1 (t) = 60 + 50 × sin(2π t/T), ν j (t) = α j + α j × sin(2π t/T), j = 2, ..., 21, ω 1,k = 20000(x + 0.001) exp(1 -500x), k = 2, ..., 11, Setting 1.2:
where α j is generated from N(30, 5 2 ).W eletthesupportsofall transfer functions be [0, 0.01], and simulate events in [0, T] with the intensity function (1) under Settings 1.1-1.2. To estimate the background intensities and transfer functions, we use cubic Bsplines with equally spaced knots. To select the numbers of Bsplines m 0 and m 1 ,wefirstperformselectionusingtheproposed procedure in Section B1 over 20 data replications. The numbers of B-splines are then fixed at the respective averages of the 20 selected values for m 0 and m 1 . It is worth noting that the estimation and selection accuracy are not overly sensitive to the number of B-splines used in the estimation (see additional results in Section B2). We have also considered larger ranges for the transfer functions and the results are very similar. We thus focus on the current setting when reporting our simulation results.
To evaluate the estimation accuracy, we report the mean squared errors. For the background intensity, it is calculated as
2 dt 1/2 and ν j (t) istheestimateofν j (t). For the transfer functions, it is calculated as MSE(ω) = 1 p j∈[p] MSE(ω j,• ),where MSE(ω j,• ) = p k=1 b 0 ( ω j,k (t)ω j,k (t)) 2 dt 1/2 and ω j,k (t) i st h ee s t i m a t eo fω j,k (t). To evaluate the selection accuracy, we report the false positive rate (FNR), the false positive rate (FPR) and the F 1 score (Forman 2003;H o ,P a r i k h , and Xing 2012), calculated as 2TP/(2TP+FP+FN), where TP is t h et ru ep o s i t i v ec o u n t ,F Pi st h ef a l s ep o s i t i v ec o u n t ,a n dF N is the false negative count. The highest F 1 score is 1 indicating perfect selection. For both NStaHawkes and StaHawkes estimators, we report the estimation and selection accuracy. For the BinGLM estimator, we only report the selection accuracy, as this method cannot be used to estimate the intensity functions. Table 1 reports the average criteria from the three methods, with standard errors in the parentheses, over 100 data replications.
TheproposedmethodNStaHawkes is seen to achieve the best performance, both in terms of the estimation accuracy and edge selection accuracy, and this holds true for different observation window length T. Moreover, it is seen that the estimation error of NStaHawkes decreases as T increases. Such an observation agrees with our theoretical result in Theorem 3.
In this simulation, we evaluate the edge selection accuracy of our proposed method. We consider two types of networks. The first type of networks are assumed to follow an Erdos-Renyi network model with edge probability p e (Erdős and Renyi 1969). In an Erdos-Renyi network model, edges are generated independently from a Bernoulli distribution with probability p e . The second type of networks follow a scale-free network model, with the degrees of nodes generated from a power-law distribution with parameter α;s u c hn e t w o r k sh a v eas k e w e d degree distributions and a larger α indicates a higher degree heterogeneity (Clauset, Shalizi, and Newman 2009). We set p = 100, p e = 0.025, 0.05 and α = 1, 2. The generated networks are shown in Figure 2. Based on the generated networks, we simulate data using the following setting.
Setting 2: ν j (t) = α j + α j × sin(2π ft/T), j = 1, ..., 100, where α j is generated independently from N(100, 5 2 ) for each node. The transfer functions are the same as in Setting 1.2. We simulate events in [0, T] with intensity function (1)w i t h f = 5a n dT = 20. To estimate the background intensities, we use cubic B-splines with equally spaced knots and to estimate the transfer functions, we use step functions with equally spaced knots, as considered in Hansen, Reynaud-Bouret, and Rivoirard (2015) and Chen et al. (2017). The numbers of basis functions m 0 and m 1 areselectedfollowingthesameprocedure as in Simulation 1.
Table 2 compares the false negative rate, false positive rate and F 1 s c o r eo ft h et h r e em e t h od so v e r1 0 0 data replications. It is seen that NStaHawkes achieves the best edge selection accuracy, in terms of F 1 scores, across all settings; StaHawkes shows a large false positive rate and this is likely due to the biased estimation of the background intensity functions; BinGLM showsalargefalsenegativerateandthisis possibly due to the loss of information in the binning approach.
We have also considered p = 200, where NStaHawkes continues to achieve a satisfactory edge selection accuracy (see Section B3).
In this section, we apply our proposed method to a neuron spike train dataset and estimate the functional connectivity network of neurons in the rat prefrontal cortex. The data were obtained from adult male Sprague-Dawley rats performing a Tmaze based delayed-alternation task of working memory (Devilbiss and Waterhouse 2004). In the experiment, the animal was t r a i n e dt on a v i g a t ed o w nt h eT -m a z ea n dc h o o s eo n eo ft w o arms (opposite to the one previously visited) for food rewards.
In each trial, the animal was released after being placed in a start box for a fixed length of delay. On a correct trial (i.e., the arm with food was chosen), the animal was rewarded and returned to the start box. On an incorrect trial (i.e., the arm without food was chosen), the animal was returned to the start box without being rewarded. In the study, the animal remained in a training period until it reached 90%-100% accuracy on 40 trials. After the training period, a recording session was performed. The spike train recording consisted of 73 neurons in an experiment of 40 trials. Each trials took about 36 sec and the total recording had 1434.22 sec. See Zhang et al. (2016) for more information about data collection and processing. We applied our proposed method to this dataset. To estimate the background intensities, we used cubic B-splines with equally spacedknotsandtoestimatethetransferfunctions,weusedstep functions with equally spaced knots, as considered in Hansen, Reynaud-Bouret, and Rivoirard (2015) and Chen et al. (2017). The numbers of basis functions m 0 and m 1 were selected using the proposed procedure in Section B1. The GIC in (9)wasused to select the tuning parameters with α T = (log p) 2 log T/2. The range of the transfer functions were set to [0, 2].W eha v ealso considered a larger range, and the results remain very similar. First, we performed the proposed test of hypothesis for each neuron to assess the if the background intensity is constant in time. Based on the p-values from the tests, 38 neurons had timevarying background intensity functions (significance level was set to 0.05). Next, we move to estimate the neuronal connectivity network using the proposed NStaHawkes. When estimating the network structure, we also considered StaHawkes and BinGLM. The tuning parameters in StaHawkes and BinGLM were selected using their recommended BIC functions, respectively. Figure 3 shows the estimated neuronal networks from the threedifferentmethods.Wecanseeallthreeestimatednetworks are sparse, with both excitatory and inhibitory relationships. However, their structures are quite different. The network estimated from our method is highly clustered and has a powerlaw degree distribution, which are two unique features of real world networks (Barabási and Albert 1999). Also interestingly, about 70% of the identified edges in our estimated network are within the right prefrontal cortex, which agrees with existing findingsthattherightprefrontalcortexishighlyrelatedtothe episodic memory retrieval (Henson, Shallice, and Dolan 1999). The biological significance of the identified edges requires further investigation.
Compared to our estimated network, StaHawkes identified a very sparse network. This difference is likely due to the bias in estimating the background intensity function from their method. BinGLM also identified a very different network structure. This network has two hub (or densely connected) nodes, namely, neurons 60 and 62, and a small clustering coefficient. We fi n dt h a tn e u r o n s6 0a n d6 2a r et h et w om o s tfr e q u e n t l yfi r e d neurons in the ensemble. Specifically, neurons 60 and 62 have