Niels Richard Hansen
email: niels.r.hansen@math.ku.dk

Patricia Reynaud-Bouret

Vincent Rivoirard
email: vincent.rivoirard@dauphine.fr

Lasso and probabilistic inequalities for multivariate point processes

Keywords: Multivariate counting process, Hawkes processes, adaptive estimation, Lasso procedure, Bernstein-type inequalities

. We observe an excellent behavior of our procedure. We rely on theoretical aspects for the essential question of tuning our methodology. Unlike adaptive Lasso of [64], our tuning procedure is proven to be robust with respect to all the parameters of the problem, revealing its potential for concrete purposes, in particular in neuroscience.

Introduction

The Lasso, proposed in [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], is a well established method that achieves sparsity of an estimated parameter vector via 1 -penalization. In this paper, we focus on using Lasso to select and estimate coefficients in the basis expansion of intensity processes for multivariate point processes.

of the distribution of ChIP-seq data along the genome [START_REF] Carstensen | Multivariate Hawkes process models of the occurrence of regulatory elements and an analysis of the pilot ENCODE regions[END_REF]. In the previous examples, the intensity of a future occurrence of a point depends on the history of all or some of the coordinates of the point processes, and it is of particular interest to estimate this dependence. This can be achieved using a parametric family of models, as in several of the papers above. Our aim is to provide a nonparametric method based on the Lasso.

The statistical properties of Lasso are particularly well understood in the context of regression with i.i.d. errors or for density estimation where a range of oracle inequalities have been established. These inequalities, now widespread in the literature, provide theoretical error bounds that hold on events with a controllable (large) probability. See for instance [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF][START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF][START_REF] Bunea | Aggregation and sparsity via 1 penalized least squares[END_REF][START_REF] Bunea | Sparse density estimation with l 1 penalties[END_REF][START_REF] Bunea | Aggregation for Gaussian regression[END_REF][START_REF] Bunea | Sparsity Oracle Inequalities for the Lasso[END_REF][START_REF] Van De Geer | High dimensional generalized linear models and the Lasso[END_REF]. We refer the reader to [START_REF] Bühlmann | Statistics for high-dimensional data[END_REF] for an excellent account on many state-of-the-art results. One main challenge in this context is to obtain as weak conditions as possible on the design -or Gram -matrix. The other important challenge is to be able to provide an 1 -penalization procedure that provides excellent performance from both theoretical and practical points of view. Standard Lasso proposed in [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] and based on deterministic constant weights constitutes a major contribution from the methodological point of view, but underestimation due to its shrinkage nature may lead to poor practical performance in some contexts. Alternative two step procedures have been suggested to overcome this drawback (see [START_REF] Meinshausen | Relaxed Lasso[END_REF][START_REF] Van De Geer | The adaptive and the thresholded Lasso for potentially misspecified models (and a lower bound for the Lasso)[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]). Zou in [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF] also discusses problems for standard Lasso to cope with variable selection and consistency simultaneously. He overcomes these problems by introducing non-constant data-driven 1 -weights based on preliminary consistent estimates.

Our contributions

In this paper we consider an 1 -penalized least squares criterion for the estimation of coefficients in the expansion of a function parameter. As in [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF][START_REF] Huang | Adaptive Lasso for sparse highdimensional regression models[END_REF][START_REF] Van De Geer | The adaptive and the thresholded Lasso for potentially misspecified models (and a lower bound for the Lasso)[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF], we consider non-constant data-driven weights. However the setup is here that of multivariate point processes and the function parameter that lives in a Hilbert space determines the point process intensities. Even in this unusual context, the least squares criterion also involves a random Gram matrix as well, and in this respect, we present a fairly standard oracle inequality with a strong condition on this Gram matrix (see Theorem 1 in Section 2).

One major contribution of this article is to provide probabilistic results that enable us to calibrate 1 -weights in the most general setup (see Theorem 2 in Section 2). This is naturally linked to sharp Bernstein type inequalities for martingales. In the literature, those kinds of inequalities generally provide upper bounds for the martingale that are deterministic and unobservable [START_REF] Shorack | Empirical processes with applications to statistics[END_REF][START_REF] Van De Geer | Exponential inequalities for martingales, with application to maximum likelihood estimation for counting processes[END_REF]. To choose data-driven weights we need observable bounds. More recently, there have been some attempts to use self-normalized processes in order to provide more flexible and random upper bounds [START_REF] Bercu | Exponential inequalities for self-normalized martingales with applications[END_REF][START_REF] De La Peña | A general class of exponential inequalities for martingales and ratios[END_REF][START_REF] De La Peña | Self-Normalized Processes[END_REF][START_REF] Dzhaparidze | On Bernstein-type inequalities for martingales[END_REF]. Nevertheless, those bounds are usually not (completely) observable when dealing with counting processes. We prove a result that goes further in this direction by providing a completely sharp random observable upper bound for the martingale in our counting process framework (see Theorem 3 in Section 3).

The second main contribution is to provide a quite theoretical and abstract framework imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

to deal with processes whose intensity is (or is well approximated by) a linear transformation of deterministic parameters to infer. This general framework also allows for different asymptotics in terms of the number of observed processes or in terms of the duration of the recording of observations, according to the setup. We focus in this paper on three main examples: the Poisson model, the Aalen multiplicative intensity model and the multivariate Hawkes process, but many other situations can be expressed in the present framework, which has the advantage of full flexibility. The first two examples have been extensively studied in the literature as we detail hereafter, but Hawkes processes are typical of situations where very little is known from a nonparametric point of view, and where fully implementable adaptive methods do not exist until the present work, to the best of our knowledge. They also constitute processes that are often used in practicein particular in neuroscience -as explained below.

It is also notable that we, in each of these three previous examples, can verify explicitly if the strong condition on the Gram matrix mentioned previously is fulfilled with probability close to 1 (see Section 4 for the Poisson and Aalen cases and Section 5 for the Hawkes case). For the multivariate Hawkes process this involves novel probabilistic inequalities. Even though the Hawkes processes have been studied extensively in the literature, see [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF][START_REF] Daley | An introduction to the theory of point processes[END_REF], very little is known about exponential inequalities and non-asymptotic tail control. Besides the univariate case [START_REF] Reynaud-Bouret | Some non asymptotic tail estimates for Hawkes processes[END_REF], no exponential inequality controlling the number of points per interval is known to us. We derive such results and other sharp controls on the convergence in the ergodic theorem to obtain control on the Gram matrix.

Finally, we carry out a simulation study in Section 6 for the most intricate process, namely the multivariate Hawkes process, with a main aim: to convince practitioners, for instance in neuroscience, that this method is indeed fully implementable and gives good results in practice. Data-driven weights for practical purposes are slight modifications of theoretical ones. These modifications essentially aim at reducing the number of tuning parameters to one. Due to non-negligible shrinkage that is unavoidable, in particular for large coefficients, we propose a two step procedure where estimation of coefficients is handled by using ordinary least squares estimation on the support preliminary determined by our Lasso methodology. Tuning issues are extensively investigated in our simulation study, and Table 1 in Section 6.3 shows that our methodology can easily and robustly be tuned by using limit values imposed by assumptions of Theorem 2. We naturally compare our procedure with adaptive Lasso of [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF] for which weights are proportional to the inverse of ordinary least squares estimates. The latter is very competitive for estimation aspects since shrinkage becomes negligible if the preliminary OLS estimates are large. But adaptive Lasso does not incorporate random fluctuations of coefficient estimators. So it is most of the time outperformed by our procedure. In particular, tuning adaptive Lasso in the Hawkes setting is a difficult task, which cannot be tackled by using standard cross-validation. Our simulation study shows that the performance of adaptive Lasso is very sensitive to the choice of the tuning parameter. Robustness with respect to tuning is another advantage of our method over adaptive Lasso. For simulations, the framework of neuronal networks is used. Our short study proves that our methodology can be used for solving concrete problems in neuroscience such as the inference of functional connectivity graphs.

Multivariate counting process

The framework introduced here and used throughout the paper aims at unifying several situations, making the reading easier. Main examples are then shortly described, illustrating the use of this setup.

We consider an M -dimensional counting process (N (m) t

) m=1,...,M , which can also be seen as a random point measure on R with marks in {1, . . . , M }, and corresponding predictable intensity processes (λ (m) t ) m=1,...,M under a probability measure P (see [START_REF] Brémaud | Point processes and queues[END_REF] or [START_REF] Daley | An introduction to the theory of point processes[END_REF] for precise definitions).

Classical models assume that the intensity λ (m) t can be written as a linear predictable transformation of a deterministic function parameter f 0 belonging to a Hilbert space H (the structure of H, and then of f 0 , will differ according to the context, as illustrated below). We denote this linear transformation by ψ(f ) = (ψ (1) (f ), ..., ψ (M ) (f )).

(1.1)

Therefore, for classical models, for any t, λ

(f 0 ).

(

The main goal in classical settings is to estimate f 0 based on observing (N (m) t

) m=1,...,M for t ∈ [0, T ]. Actually, we do not require in Theorems 1 and 2 that (1.2) holds. Our aim is mainly to furnish an estimate of the best linear approximation ψ (m) t (f 0 ) of the underlying intensity λ (m) t . Let us illustrate the general setup with three main examples: First, the case with i.i.d. observations of an inhomogeneous Poisson process on [0, 1] and unknown intensity, second, the well known Aalen multiplicative intensity model and third, the central example of the multivariate Hawkes process. For the first two models, asymptotics is with respect to M (T is fixed). For the third one, M is fixed and asymptotics is with respect to T .

The Poisson model

Let us start with a very simple example which will be somehow a toy problem here compared to the other examples. In this example we take T = 1 and assume that we observe M i.i.d. Poisson processes on [0, 1] with common intensity f 0 : [0, 1] -→ R + . Asymptotic properties are obtained when M tends to infinity. In this case, the intensity λ (m) of the m'th process N (m) is f 0 , which does not depend on m: Therefore for any m ∈ {1, . . . , M } and any t, we set This framework has already been extensively studied from an adaptive point of view: see for instance [START_REF] Reynaud-Bouret | Adaptive estimation of the intensity of inhomogeneous Poisson processes via concentration inequalities[END_REF][START_REF] Willett | Multiscale Poisson Intensity and Density Estimation[END_REF] for model selection methods, [START_REF] Reynaud-Bouret | Near optimal thresholding estimation of a Poisson intensity on the real line[END_REF] for wavelet thresholding estimation or [START_REF] Reynaud-Bouret | Spike trains as (in)homogeneous Poisson processes or Hawkes processes: nonparametric adaptive estimation and goodness-of-fit tests[END_REF] for kernel estimates. In this context, our present general result matches with existing minimax adaptation results where asymptotics is with respect to M .

ψ (m) t (f 0 ) := f 0 (t), and H = L 2 ([0, 1]) is equipped with the classical norm defined by ||f || = 1 0 f 2 (t)dt

The Aalen multiplicative intensity model

This is one of the most popular counting processes because of its adaptivity to various situations (from Markov models to censored lifetime models) and its various applications to biomedical data (see [START_REF] Andersen | Statistical models based on counting processes[END_REF]). Given X a Hilbert space, we consider f 0 : [0, T ] × X -→ R + , and we set for any t ∈ R,

λ (m) t = ψ (m) t (f 0 ) := f 0 (t, X (m) )Y (m) t ,
where Y (m) is an observable predictable process and X (m) represents covariates. In this case, H = L 2 ([0, T ] × X ). Our goal is to estimate f 0 and not to select covariates. So, to fix ideas one sets X = [0, 1] and T = 1. Hence H can be identified with L 2 ([0, 1] 2 ). For right-censored data, f 0 usually represents the hazard rate. The presence of covariates in this pure nonparametric model is the classical generalization of the semi-parametric model proposed by Cox (see [START_REF] Letue | Modèle de Cox : Estimation par sélection de modèle et modèle de chocs bivarié[END_REF] for instance). Note that the Poisson model is a special case of the Aalen model.

The classical framework consists in assuming that (X (m) , Y (m) , N (m) ) m=1,...,M is an i.i.d. M -sample and as for the Poisson model, it is natural to investigate asymptotic properties when M → +∞. If there are no covariates, several adaptive approaches already exist: See [START_REF] Brunel | Penalized contrast estimation of density and hazard rate with censored data[END_REF][START_REF] Brunel | Adaptive estimation of hazard rate with censored data[END_REF][START_REF] Reynaud-Bouret | Penalized projection estimators of the Aalen multiplicative intensity[END_REF] for various penalized least-squares contrasts and [START_REF] Chagny | Adaptive warped kernel estimators[END_REF] for kernel methods in special cases of censoring. In the presence of covariates, one can mention [START_REF] Aalen | A model for nonparametric regression analysis of counting processes[END_REF][START_REF] Andersen | Statistical models based on counting processes[END_REF] for a parametric approach, [START_REF] Comte | Adaptive estimation of the conditional intensity of marker-dependent counting processes[END_REF][START_REF] Letue | Modèle de Cox : Estimation par sélection de modèle et modèle de chocs bivarié[END_REF] for a model selection approach and [START_REF] Gaïffas | High-dimensional additive hazard models and the Lasso[END_REF] for a Lasso approach. Let us also cite [START_REF] Bunea | Covariate selection for semiparametric hazard function regression models[END_REF] where covariate selection via penalized MLE has been studied. Once again, our present general result matches with existing oracle results. In [START_REF] Chagny | Adaptive warped kernel estimators[END_REF], exponential control of random fluctuations leading to adaptive results are derived without using the martingale theory. In more general frameworks (as in [START_REF] Comte | Adaptive estimation of the conditional intensity of marker-dependent counting processes[END_REF] for instance), martingales are required and this even when i.i.d. processes are involved.

Hawkes processes

Hawkes processes are the point processes equivalent to autoregressive models. In seismology, Hawkes processes can model earthquakes and their aftershocks [START_REF] Vere-Jones | Some examples of statistical estimation applied to earthquake data I: cyclic Poisson and self-exciting models[END_REF]. More recently they have been used to model favored or avoided distances between occurrences of motifs [START_REF] Gusto | FADO: a statistical method to detect favored or avoided distances between occurrences of motifs using the Hawkes' model[END_REF] or Transcription Regulatory Elements [START_REF] Carstensen | Multivariate Hawkes process models of the occurrence of regulatory elements and an analysis of the pilot ENCODE regions[END_REF] on the DNA. We can also mention the use of Hawkes processes as models of social interactions [START_REF] Mitchell | Hawkes process as a model of social interactions: a view on video dynamics[END_REF] or financial phenomena [START_REF] Bacry | Some limit theorems for Hawkes processes and application to financial statistics[END_REF].

In the univariate setting, with M = 1, the intensity of a nonlinear Hawkes process (N t ) t>0 is given by

λ t = φ t- -∞ h(t -u)dN u ,
imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: [START_REF] Chagny | Adaptive warped kernel estimators[END_REF] where φ : R → R + and h : R + → R (see [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF]). A particular case is Hawkes's self exciting point process, for which h is nonnegative and φ(x) = ν + x where ν > 0 (see [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF][START_REF] Daley | An introduction to the theory of point processes[END_REF][START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF]). For instance, for seismological purposes, ν represents the spontaneous occurrences of real original earthquakes. The function h models self-interaction: after a shock at time u, we observe an aftershock at time t with large probability if h(t -u) is large.

These notions can be easily extended to the multivariate setting and in this case the intensity of N (m) takes the form:

λ (m) t = φ (m) M =1 t- -∞ h (m) (t -u)dN ( ) (u) .
Theorem 7 of [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF] gives conditions on the functions φ (m) (namely Lipschitz properties) and on the functions h (m) to obtain existence and uniqueness of a stationary version of the associated process. Throughout this paper, we assume that for any m ∈ {1, . . . , M },

φ (m) (x) = (ν (m) + x) + ,
where ν (m) > 0 and (•) + denotes the positive part. Note that in [START_REF] Carstensen | Multivariate Hawkes process models of the occurrence of regulatory elements and an analysis of the pilot ENCODE regions[END_REF][START_REF] Chornoboy | Maximum likelihood identification of neural point process systems[END_REF], the case φ (m) (x) = exp(ν (m) +x) was studied. However, Lipschitz properties required in [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF] are not satisfied in this case. By introducing, as previously, the linear predictable transformation ψ(f ) = (ψ (1) (f ), ..., ψ (M ) (f )) with for any m and any t ψ

(m) t (f 0 ) := ν (m) + M =1 t- -∞ h (m) (t -u)dN ( ) (u), (1.3)
with f 0 = (ν (m) , (h (m) ) =1,...,M ) m=1,...,M , we have λ

(m) t = (ψ (m) t (f 0 )) + .
Note that the upper integration limits in (1.3) are t-, that is, the integrations are all over the open interval (-∞, t). This assures predictability of the intensity disregarding the values of h (m) (0). Alternatively, it can be assumed throughout that h (m) (0) = 0, in which case the integrals in (1.3) can be over (-∞, t] without compromising predictability. The parameter ν (m) is called the spontaneous rate, whereas the function h (m) is called the interaction function of N ( ) on N (m) . The goal is to estimate f 0 by using Lasso estimates. In the sequel, we will assume that the support of h (m) is bounded. By rescaling we can then assume that the support is in [0, 1], and we will do so throughout. Note that in this case we will need to observe the process on [-1, T ] in order to compute ψ

(m) t (f 0 ) for t ∈ [0, T ].
The Hilbert space H associated with this setup is

H = (R × L 2 ([0, 1]) M ) M = f = (µ (m) , (g (m) ) =1,...,M ) m=1,...,M : g (m) with support in [0, 1] and ||f || 2 = m (µ (m) ) 2 + m 1 0 g (m) (t)
λ (m) t = ψ (m) t (f 0 ).
The multivariate point process associated with this setup is called the multivariate Hawkes self exciting point process (see [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF]). In this example, M is fixed and asymptotic properties are obtained when T tends to infinity.

To the best of our knowledge, nonparametric estimation for Hawkes models has only been proposed in [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF] in the univariate setting where the method is based on 0penalization of the least-squares contrast. However, due to 0 -penalization, the criterion is not convex and the computational cost, in particular for the memory storage of all the potential estimators, is huge. Therefore, this method has never been adapted to the multivariate setting. Moreover, the penalty term in this method is not data-driven and ad-hoc tuning procedures have been used for simulations. This motivates the present work and the use of a convex Lasso criterion combined with data-driven weights, to provide a fully implementable and theoretically valid data-driven method, even in the multivariate case.

Applications in neuroscience.

Hawkes processes can naturally be applied to model neuronal activity. Extracellular action potentials can be recorded by electrodes and the recorded data for the neuron m can be seen as a point process, each point corresponding to the peak of one action potential of this neuron (see [START_REF] Brette | Handbook of Neural Activity Measurement[END_REF] for instance for precise definitions). When M neurons are simultaneously recorded, one can assume that we are faced with a realization of N = (N (m) ) m=1,...,M modeled by a multivariate Hawkes process. We then assume that the intensity associated with the activity of the neuron m is given by λ

(m) t = (ψ (m) t (f 0 )) + , where ψ (m) t (f 0 ) is given in (1.3). At any occurrence u < t of N ( ) , ψ (m) t (f 0 ) increases (excitation) or decreases (inhibition) according to the sign of h (m) (t -u).
Modeling inhibition is essential from the neurobiological point of view. So, we cannot assume that all interaction functions are nonnegative, and we cannot omit the positive part. More details are given in Section 6.

In neuroscience, Hawkes processes combined with maximum likelihood estimation have been used in the seminal paper [START_REF] Chornoboy | Maximum likelihood identification of neural point process systems[END_REF], but the application of the method requires a too huge number of observations for realistic practical purposes. Models based on Hawkes processes have nevertheless been recently discussed in neuroscience, since they constitute one of the few simple models able to produce dependence graphs between neurons, that may be interpreted in neuroscience as functional connectivity graphs [START_REF] Pernice | How structure determines correlations in neuronal networks[END_REF][START_REF] Pernice | Recurrent interactions in spiking networks with arbitrary topology[END_REF]. However, many nonparametric statistical questions arise that are not solved yet in order to furnish a fully applicable tool for real data analysis [START_REF] Krumin | Correlation-based analysis and generation of multiple spike trains using Hawkes models with an exogenous input[END_REF]. We think that the Lasso-based methodology presented in this paper may furnish the first robust tool in this direction.

Notation and overview of the paper

Some notation from the general theory of stochastic integration is useful to simplify the otherwise quite heavy notation. If H = (H (1) , ..., H (M ) ) is a multivariate process with locally bounded coordinates, say, and X = (X (1) , ..., X (M ) ) is a multivariate semimartingale, we define the real valued process H • X by

H • X t := M m=1 t 0 H (m) s dX (m) s .
Given F : R -→ R we use F(H) to denote the coordinatewise application of F, that is F(H) = (F(H (1) ), ..., F(H (M ) )). In particular,

F(H) • X t = M m=1 t 0 F(H (m) s )dX (m) s .
We also define the following scalar product on the space of multivariate processes. For any multivariate processes H = (H (1) , ..., H (M ) ) and K = (K (1) , ..., K (M ) ), we set

< H, K > proc := M m=1 T 0 H (m) s K (m)
s ds, the corresponding norm being denoted ||H|| proc . Since ψ introduced in (1.1) is linear, the Hilbert space H inherits a bilinear form from the previous scalar product, that we denote, for all f, g in H,

< f, g > T :=< ψ(f ), ψ(g) > proc = M m=1 T 0 ψ (m) s (f )ψ (m) s (g)ds,
and the corresponding quadratic form is denoted ||f || 2 T . The compensator Λ = (Λ (m) ) m=1,...,M of N = (N (m) ) m=1,...,M is finally defined for all t by

Λ (m) t = t 0 λ (m)
s ds.

Section 2 gives our main oracle inequality and the choice of the 1 -weights in the general framework of counting processes. Section 3 provides the fundamental Bernsteintype inequality. Section 4 details the meaning of the oracle inequality in the Poisson and Aalen setups. The probabilistic results needed for the Hawkes processes as well as the interpretation of the oracle inequality in this framework is done in Section 5. Simulations on multivariate Hawkes processes are performed in Section 6. The last Section is dedicated to the proofs of our results.

imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

Lasso estimate and oracle inequality

We wish to estimate the true underlying intensity so our main goal consists in estimating the parameter f 0 . For this purpose we assume we are given Φ a dictionary of functions (whose cardinality is denoted |Φ|) and we define f a as a linear combination of the functions of Φ, that is,

f a := ϕ∈Φ a ϕ ϕ,
where a = (a ϕ ) ϕ∈Φ belongs to R Φ . Then, since ψ is linear, we get

ψ(f a ) = ϕ∈Φ a ϕ ψ(ϕ).
We use the following least-squares contrast C defined on H by

C(f ) := -2 ψ(f ) • N T + ||f || 2 T , ∀ f ∈ H. (2.1)
This contrast, or some variations of C, have already been used in particular setups (see for instance [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF] or [START_REF] Gaïffas | High-dimensional additive hazard models and the Lasso[END_REF]). The main heuristic justification lies in following arguments. Since ψ(f ) is a predictable process, the compensator at time T of C(f ) is given by

C(f ) := -2ψ(f ) • Λ T + ||f || 2 T ,
which can also be written as C(f ) = -2 < ψ(f ), λ > proc +||ψ(f )|| 2 proc . Note that C is minimum when ||ψ(f ) -λ|| proc is minimum. If λ = ψ(f 0 ) and if ||.|| T is a norm on the Hilbert space H then the unique minimizer of C is f 0 . Therefore, to get the best linear approximation of λ of the form ψ(f ), it is natural to look at minimizers of C(f ). Restricting to linear combinations of functions of Φ, since ψ is linear, we obtain

C(f a ) = -2a b + a Ga
where a denotes the transpose of the vector a and for ϕ 1 , ϕ 2 ∈ Φ,

b ϕ1 = ψ(ϕ 1 ) • N T , G ϕ1,ϕ2 =< ϕ 1 , ϕ 2 > T .
(2.2) Note that both the vector b of dimension |Φ| and the Gram matrix G of dimensions |Φ| × |Φ| are random but nevertheless observable.

To estimate a we minimize the contrast, C(f a ), subject to an 1 -penalization on the a-vector. That is, we introduce the following 1 -penalized estimator

â ∈ argmin a∈R Φ {-2a b + a Ga + 2d |a|} (2.3)
where |a| = (|a ϕ |) ϕ∈Φ and d ∈ R Φ + . With a good choice of d the solution of (2.3) will achieve both sparsity and good statistical properties. Finally, we let f = f â denote the Lasso estimate associated with â.

Our first result establishes theoretical properties of f by using the classical oracle approach. More precisely, we establish a bound on the risk of f if some conditions are true. This is a non-probabilistic result that only relies on the definition of â by (2.3). In the next section we will deal with the probabilistic aspect, which is to prove that the conditions are fulfilled with large probability.

Theorem 1. Let c > 0. If G cI (2.4
)

and if for all ϕ ∈ Φ |b ϕ -bϕ | ≤ d ϕ , (2.5
)

where bϕ = ψ(ϕ) • Λ T ,
then there exists an absolute constant C, independent of c, such that

||ψ( f ) -λ|| 2 proc ≤ C inf a∈R Φ    ||λ -ψ(f a )|| 2 proc + c -1 ϕ∈S(a) d 2 ϕ    , (2.6)
where S(a) is the support of a. If λ = ψ(f 0 ), the oracle inequality (2.6) can also be rewritten as

|| f -f 0 || 2 T ≤ C inf a∈R Φ    ||f 0 -f a || 2 T + c -1 ϕ∈S(a) d 2 ϕ    , (2.7)
The proof of Theorem 1 is given in Section 7.1. Note that Assumption (2.4) ensures that G is invertible and then coordinates of â are finite almost surely. Assumption (2.4) also ensures that ||f || T is a real norm on f at least when f is a linear combination of the functions of Φ.

Two terms are involved on the right hand sides of (2.6) and (2.7). The first one is an approximation term and the second one can be viewed as a variance term providing a control of the random fluctuations of the b ϕ 's around the bϕ 's. Note that b ϕ -bϕ = ψ(ϕ)•(N -Λ) T is a martingale (see also the comments after Theorem 2 for more details). The approximation term can be small but the price to pay may be a large support of a, leading to large values for the second term. Conversely, a sparse a leads to a small second term. But in this case the approximation term is potentially larger. Note that if the function f 0 can be approximated by a sparse linear combination of the functions of Φ, then we obtain a sharp control of || f -f 0 || 2 T . In particular, if f 0 can be decomposed on the dictionary, so we can write f 0 = f a0 for some a 0 ∈ R Φ , then (2.7) gives

|| f -f 0 || 2 T ≤ Cc -1 ϕ∈S(a0) d 2 ϕ .
In this case, the right hand side can be viewed as the sum of the estimation errors made by estimating the components of a 0 .

Such oracle inequalities are now classical in the huge literature of Lasso procedures. See for instance [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF][START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF][START_REF] Bunea | Aggregation and sparsity via 1 penalized least squares[END_REF][START_REF] Bunea | Sparse density estimation with l 1 penalties[END_REF][START_REF] Bunea | Aggregation for Gaussian regression[END_REF][START_REF] Bunea | Sparsity Oracle Inequalities for the Lasso[END_REF][START_REF] Koltchinskii | Nuclear norm penalization and optimal rates for noisy low rank matrix completion[END_REF][START_REF] Van De Geer | High dimensional generalized linear models and the Lasso[END_REF], who established oracle inequalities in the same spirit as in Theorem 1. We bring out the paper [START_REF] Candes | The Dantzig selector: statistical estimation when p is much larger than n[END_REF], which gives technical and heuristic arguments for justifying optimality of such oracle inequalities (see Section 1.3 of [START_REF] Candes | The Dantzig selector: statistical estimation when p is much larger than n[END_REF]). Most of these papers deal with independent data.

In the sequel, we prove that Assumption (2.4) is satisfied with large probability by using the same approach as [START_REF] Rudelson | Smallest singular value of a random rectangular matrix[END_REF][START_REF] Rudelson | Non-asymptotic theory of random matrices: extreme singular values[END_REF] and to a lesser extent as Section 2.1 of [START_REF] Candes | The Dantzig selector: statistical estimation when p is much larger than n[END_REF] or [START_REF] Rudelson | On sparse reconstruction from Fourier and Gaussian measurements[END_REF]. Section 5 is in particular mainly devoted to show that (2.4) holds with large probability for the multivariate Hawkes processes.

For Theorem 1 to be of interest, the condition on the martingale, condition (2.5), needs to hold with large probability as well. From this control, we deduce convenient datadriven 1 -weights that are the key parameters of our estimate. Note that our estimation procedure does not depend on the value of c in (2.4). So knowing the latter is not necessary for implementing our procedure. Therefore, one of the main contributions of the paper is to provide new sharp concentration inequalities that are satisfied by multivariate point processes. This is the main goal of Theorem 3 in Section 3 where we establish Bernstein type inequalities for martingales. We apply it to the control of (2.5). This allows us to derive the following result, which specifies the choice of the d ϕ 's needed to obtain the oracle inequality with large probability.

Ω c = {G cI} .
Let x and ε be strictly positive constants and define for all ϕ ∈ Φ,

d ϕ = 2(1 + ε) V µ ϕ x + B ϕ x 3 , (2.8
)

with V µ ϕ = µ µ -φ(µ) (ψ(ϕ)) 2 • N T + B 2 ϕ x µ -φ(µ)
for a real number µ such that µ > φ(µ), where φ(u) = exp(u) -u -1. Let us consider the Lasso estimator f of f 0 defined in Section 2. Then, with probability larger than

1 -4 ϕ∈Φ   log 1 + µVϕ B 2 ϕ x log(1 + ε) + 1   e -x -P(Ω c V,B ) -P(Ω c c ), inequality (2.7) is satisfied, i.e. ||ψ( f ) -λ|| 2 proc ≤ C inf a∈R Φ    ||λ -ψ(f a )|| 2 proc + c -1 ϕ∈S(a) d 2 ϕ    . If moreover λ = ψ(f 0 ), then || f -f 0 || 2 T ≤ C inf a∈R Φ    f 0 -f a 2 T + c -1 ϕ∈S(a) d 2 ϕ    ,
where C is a constant independent of c, Φ, T and M.

The first oracle inequality gives a control of the difference between the true intensity and ψ( f ). The equality λ = ψ(f 0 ) is not required and we can apply this result, for instance, with λ = (ψ(f 0 )) + .

Of course, the smaller the d ϕ 's the better the oracle inequality. Therefore, when x increases, the probability bound and the d ϕ 's increase and we have to realize a compromise to obtain a meaningful oracle inequality on an event with large probability. The choice of x is deeply discussed below, in Sections 4 and 5 for theoretical purposes and in Section 6 for practical purposes.

Let us first discuss more deeply the definition of d ϕ (derived from subsequent Theorem 3) which seems intricate. Up to a constant depending on the choice of µ and ε, d ϕ is of same order as max

x(ψ(ϕ)) 2 • N T , B ϕ x . To give more insight on the values of d ϕ , let us consider the very special case where for any m ∈ {1, . . . , M } for any s, ψ (m) s (ϕ) = c m 1 {s∈Am} , where c m is a positive constant and A m a compact set included into [0, T ]. In this case, by naturally choosing B ϕ = max 1≤m≤M c m , we have:

x(ψ(ϕ)) 2 • N T ≥ B ϕ x ⇐⇒ M m=1 c 2 m N (m) Am ≥ x max 1≤m≤M c 2 m ,
where

N (m)
Am represents the number of points of N (m) falling in A m . For more general vector functions ψ(ϕ), the term x(ψ(ϕ)) 2 • N T will dominate B ϕ x if the number of points of the process lying where ψ(ϕ) is large, is significant. In this case, the leading term in d ϕ is expected to be the quadratic term 2(1 + ε) µ µ-φ(µ) x(ψ(ϕ)) 2 • N T and the linear terms in x can be viewed as residual terms. Furthermore, note that when µ tends to 0,

µ µ -φ(µ) = 1 + µ 2 + o(µ), x µ -φ(µ) ∼ x µ → +∞ since x > 0.
So, if µ and ε tend to 0, the quadratic term tends to 2x(ψ(ϕ)) 2 • N T , but the price to pay is the explosion of the linear term in x. In any case, it is possible to make the quadratic term as close to 2x(ψ(ϕ)) 2 • N T as desired. Basically, this term cannot be improved (see comments after Theorem 3 for probabilistic arguments).

imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

Let us now discuss the choice of x. In more classical contexts such as density estimation based on an n-sample, the choice x = γ log(n) plugged in the parameters analog to the d ϕ 's is convenient, since it both ensures a small probability bound and a meaningful order of magnitude for the oracle bound (see [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF] for instance). See also Sections 4 and 5 for similar evaluations in our setup. But it has also been further established that the choice γ = 1 is the best. Indeed if the components of d are chosen smaller than the analog of 2x(ψ(ϕ)) 2 • N T in the density framework, then the resulting estimator is definitely bad from the theoretical point of view, but simulations also show that, to some extent, if the components of d are larger than the analog of 2x(ψ(ϕ)) 2 • N T , then the estimator deteriorates too. A similar result is out of reach in our setting, but similar conclusions may remain valid here since density estimation often provides some clues about what happens for more intricate heteroscedastic models. In particular, the main heuristic justifying the optimality of this tuning result in the density setting is that the quadratic term (and in particular the constant √ 2) corresponds to the rate of the central limit theorem and in this sense, it provides the "best approximation" for the fluctuations. For further discussion, see the simulation study in Section 6.

Finally, it remains to control P(Ω V,B ) and P(Ω c ). These are the goals of Section 4 for Poisson and Aalen models and Section 5 for multivariate Hawkes processes.

Bernstein type inequalities for multivariate point processes

We establish a Bernstein type concentration inequality based on boundedness assumptions. This result, which has an interest per se from the probabilistic point of view, is the key result to derive the convenient values for the vector d in Theorem 2 and so is capital from the statistical perspective. with respect to some given filtration. Let B > 0. Let H = (H (m) ) m=1,...,M be a multivariate predictable process such that for all ξ ∈ (0, 3), for all t,

exp(ξH/B) • Λ t < ∞ a.s. and exp(ξH 2 /B 2 ) • Λ t < ∞ a.s. (3.1)
Let us consider the martingale defined for all t ≥ 0 by

M t = H • (N -Λ) t .
Let v > w and x be positive constants and let τ be a bounded stopping time. Let us define

V µ = µ µ -φ(µ) H 2 • N τ + B 2 x µ -φ(µ)
for a real number µ ∈ (0, 3) such that µ > φ(µ), where φ(u) = exp(u) -u -1. Then, for any ε > 0,

P Mτ ≥ 2(1 + ε) V µ x + Bx 3 and w ≤ V µ ≤ v and sup m,t≤τ |H (m) t | ≤ B ≤ 2 log(v/w) log(1 + ε) + 1 e -x . (3.2)
This result is based on the exponential martingale for counting processes, which has been used for a long time in the context of the counting process theory. See for instance [START_REF] Brémaud | Point processes and queues[END_REF], [START_REF] Shorack | Empirical processes with applications to statistics[END_REF] or [START_REF] Van De Geer | Exponential inequalities for martingales, with application to maximum likelihood estimation for counting processes[END_REF]. This basically gives a concentration inequality taking the following form (the result is stated here in its univariate form for comparison purposes): for any x > 0,

P M τ ≥ 2ρx + Bx 3 and τ 0 H 2 s λ(s)ds ≤ ρ and sup s∈[0,τ ] |H s | ≤ B ≤ e -x . (3.3)
In (3.3), ρ is a deterministic upper bound of v = τ 0 H 2 s λ(s)
ds, the bracket of the martingale, and consequently the martingale equivalent of the variance term. Moreover B is a deterministic upper bound of sup s∈[0,τ ] |H s |. The leading term for moderate values of x and τ large enough is consequently √ 2ρx. The central Limit Theorem for martingales states that, under some assumptions, a sequence of martingales (M n ) n with respective brackets (v n ) n tending to a deterministic value v, once correctly normalized, tends to a Gaussian process with bracket v. Therefore, a term of the form √ 2vx in the upper bound is not improvable, in particular the constant √ 2. However the replacement of the limit v by a deterministic upper bound ρ constitutes a loss. In this sense, Theorem 3 improves the bound and consists of plugging in the unbiased estimate v = τ 0 H 2 s dN s instead of a non sharp deterministic upper bound of v. Note that we are not able to obtain exactly the term √ 2 but any value strictly larger than √ 2, as close as we want to √ 2 up to some additive terms depending on B that are negligible for moderate values of x.

The proof is based on a peeling argument that was first introduced in [START_REF] Liptser | Deviation probability bound for martingales with applications to statistical estimation[END_REF] for Gaussian processes and is given in Section 7.3.

Note that there exist also inequalities that seem nicer than (3.3) which constitutes the basic brick for our purpose. For instance, in [START_REF] Dzhaparidze | On Bernstein-type inequalities for martingales[END_REF], it is established that for any deterministic positive real number θ, for any x > 0,

P M τ ≥ √ 2θx and τ 0 H 2 s dΛ s + τ 0 H 2 s dN s ≤ θ ≤ e -x . (3.4)
At first sight, this seems better than Theorem 3 because no linear term depending on B appears, but if we wish to use the estimate 2 τ 0 H 2 s dN s instead of θ in the inequality, we have to bound |H s | by some B in any case. Moreover, by doing so, the quadratic term will be of order √ 4vx which is worse than the term √ 2vx derived in Theorem 3, even if this constant √ 2 can only be reached asymptotically in our case.

There exists a better result if the martingale M t is for instance conditionally symmetric (see [START_REF] De La Peña | A general class of exponential inequalities for martingales and ratios[END_REF][START_REF] Bercu | Exponential inequalities for self-normalized martingales with applications[END_REF][START_REF] Dzhaparidze | On Bernstein-type inequalities for martingales[END_REF]): for any x > 0,

P M τ ≥ √ 2κx and τ 0 H 2 s dN s ≤ κ ≤ e -x , (3.5)
which seems close to the ideal inequality. But there are actually two major problems with this inequality. First, one needs to assume that the martingale is conditionally symmetric, which cannot be the case in our situation for general counting processes and general dictionaries. Second, it depends on the deterministic upper bound κ instead of v. To replace κ by v and then to apply peeling arguments as in the proof of Theorem 3, we need to assume the existence of a positive constant w such that v ≥ w. But if the process is empty, then v = 0, so we cannot generally find such a positive lower bound, whereas in our theorem, we can always take w = B 2 x µ-φ(µ) as a lower bound for V µ . Finally, note that in Proposition 6 (see Section 7.3), we also derive a similar bound where V µ is replaced by τ 0 H 2 s dΛ s . Basically, it means that the same type of results hold for the quadratic characteristic instead of the quadratic variation. Though this result is of little use here, since the quadratic characteristic is not observable, we think that it may be of interest for readers investigating self-normalized results as in [START_REF] De La Peña | Self-Normalized Processes[END_REF].

Applications to the Poisson and Aalen models

We now apply Theorem 2 to the Poisson and Aalen models. The case of the multivariate Hawkes process, which is much more intricate, will be the subject of the next section.

The Poisson model

Let us recall that in this case, we observe M i.i.d. Poisson processes with intensity f 0 supported by [0, 1] (with M ≥ 2) and that the norm is given by ||f || 2 = 1 0 f 2 (x)dx. We assume that Φ is an orthonormal system for ||.||. In this case,

||.|| 2 1 = M ||.|| 2 and G = M I,
where I is the identity matrix. One applies Theorem 2 with c = M (so P(Ω c c ) = 0) and

B ϕ = ||ϕ|| ∞ , V ϕ = ||ϕ|| 2 ∞ (1 + δ)M m 1 ,
for δ > 0 and m 1 = 1 0 f 0 (t)dt. Note that here T = 1 and therefore

N (m) T = N (m) 1
is the total number of observed points for the mth process. Using is the Poisson distribution with parameter M m 1 , Cramer-Chernov arguments give:

ψ(ϕ) 2 • N 1 ≤ ||ϕ|| 2
P(Ω c V,B ) ≤ P M m=1 N (m) 1 > (1 + δ)M m 1 ≤ exp (-{(1 + δ) ln(1 + δ) -δ}M m 1 ) .
For α > 0, by choosing x = α log(M ), we finally obtain the following corollary derived from Theorem 2.

Corollary 1. With probability larger than

1 -C 1 |Φ| log(M ) M α -e -C2M
, where C 1 is a constant depending on µ, ε, α, δ and m 1 and C 2 is a constant depending on δ and m 1 , we have:

|| f -f0|| 2 ≤ C inf a∈R Φ    f0 -fa 2 + 1 M 2 ϕ∈S(a) log(M ) M m=1 1 0 ϕ 2 (x)dN (m) x + log 2 (M )||ϕ|| 2 ∞    ,
where C is a constant depending on µ, ε, α, δ and m 1 .

To shed some lights on this result, consider an asymptotic perspective by assuming that M is large. Assume also, for sake of simplicity, that f 0 is bounded from below on [0, 1]. If the dictionary Φ (whose size may depend on M ) satisfies

max ϕ∈Φ ||ϕ|| ∞ = o M log M , then, since, almost surely, 1 M M m=1 1 0 ϕ 2 (x)dN (m) x M →∞ -→ 1 0 ϕ 2 (x)f 0 (x)dx, almost surely, 1 M 2 ϕ∈S(a) log(M ) M m=1 1 0 ϕ 2 (x)dN (m) x + log 2 (M )||ϕ|| 2 ∞ = log M M ϕ∈S(a) 1 0 ϕ 2 (x)f 0 (x)dx × (1 + o(1)).
The right hand term corresponds, up to the logarithmic term, to the sum of variance terms when estimating

1 0 ϕ(x)f 0 (x)dx with 1 M M m=1 1 0 ϕ(x)dN (m) x
for ϕ ∈ S(a). This means that the estimator adaptively achieves the best trade-off between a bias term and a variance term. The logarithmic term is the price to pay for adaptation. Furthermore, when M → +∞, the inequality of Corollary 1 holds with probability that goes to 1 at a polynomial rate. We refer the reader to [START_REF] Reynaud-Bouret | Near optimal thresholding estimation of a Poisson intensity on the real line[END_REF] for a deep discussion on optimality of such results.

The Aalen model

Results similar to those presented in this paragraph can be found in [START_REF] Gaïffas | High-dimensional additive hazard models and the Lasso[END_REF] under restricted eigenvalues conditions instead of (2.4). Recall that we observe an M -sample (X (m) , Y (m) , N (m) ) m=1,...,M , with

Y (m) = (Y (m) t ) t∈[0,1] and N (m) = (N (m) t ) t∈[0,1] (with M ≥ 2). We assume that X (m) ∈ [0, 1] and that the intensity of N (m) t is f 0 (t, X (m) )Y (m) t . We set for any f , ||f || 2 e := E 1 0 f 2 (t, X (1) ) Y (1) t 2 dt .
We assume that Φ is an orthonormal system for || • || 2 , the classical norm on L 2 ([0, 1] 2 ), and we assume that there exists a positive constant r such that

∀ f ∈ L 2 ([0, 1] 2 ), ||f || e ≥ r||f || 2 , (4.1)
so that || • || e is a norm. If we assume, for instance, that the density of the X (m) 's is lower bounded by a positive constant c 0 and there exists c 1 > 0 such that for any t,

E[(Y (1)
t ) 2 |X (1) ] ≥ c 1 then (4.1) holds with r 2 = c 0 c 1 . The empirical version of ||f || e , denoted ||f || emp , is defined by ||f || 2 emp := 1 M ||f || 2 T = 1 M M m=1 1 0 f 2 (t, X (m) ) Y (m) t 2 dt.
Unlike the Poisson model, since the intensity depends on covariates X (m) 's and variables Y (m) 's, the control of P(Ω c c ) is much more cumbersome for the Aalen case, even if it is less intricate than for Hawkes processes (see Section 5). We have the following result proved in Section 7.5.1.

Proposition 1. We assume that (4.1) is satisfied, the density of the covariates X (m) is bounded by D and

sup t∈[0,1] max m∈{1,...,M } Y (m) t ≤ 1 almost surely. (4.2)
We consider an orthonormal dictionary Φ of functions of L 2 ([0, 1] 2 ) that may depend on M , and we let r Φ denote the spectral radius of the matrix H whose components are

H ϕ,ϕ = |ϕ(t, x)||ϕ (t, x)|dtdx. Then, if max ϕ∈Φ ||ϕ|| 2 ∞ × r Φ |Φ| × log M M → 0, (4.3)
when M → +∞ then, for any β > 0, there exists C 1 > 0 depending on β, D and f 0 such that with c = C 1 M , we have

P(Ω c c ) = O(|Φ| 2 M -β ).
Assumption (4.2) is usually satisfied in most of the practical examples where Aalen models are involved. See [START_REF] Andersen | Statistical models based on counting processes[END_REF] for explicit examples and see for instance [START_REF] Grégoire | Least squares cross-validation for counting process intensities[END_REF][START_REF] Reynaud-Bouret | Penalized projection estimators of the Aalen multiplicative intensity[END_REF] for other articles where this assumption is made. In the sequel, we also assume that there exists a positive constant R such that

max m∈{1,...,M } N (m) 1 ≤ R a.s. (4.4)
This assumption, considered by [START_REF] Reynaud-Bouret | Penalized projection estimators of the Aalen multiplicative intensity[END_REF], is obviously satisfied in survival analysis where there is at most one death per individual. It could have been relaxed in our setting, by considering exponential moments assumptions, to include Markov cases for instance. However this much simpler assumption allows us to avoid tedious and unnecessary technical aspects since we only wish to illustrate our results in a simple framework. Under (4.2) and (4.4), almost surely,

ψ(ϕ) 2 •N T = M m=1 1 0 [Y (m) t ] 2 ϕ 2 (t, X (m) )dN (m) t ≤ M m=1 1 0 ϕ 2 (t, X (m) )dN (m) t ≤ M R||ϕ|| 2 ∞ .
So, we apply Theorem 2 with

B ϕ = ||ϕ|| ∞ , V ϕ = M R||ϕ|| 2 ∞ (so P(Ω V,B ) = 1
) and x = α log(M ) for α > 0. We finally obtain the following corollary.

Corollary 2. Assume that (4.2) and (4.4) are satisfied. With probability larger than

1 -C 1 |Φ| log(M ) M α -P(Ω c c )
, where C 1 is a constant depending on µ, ε, α and R, we have:

|| f -f0|| 2 emp ≤ C inf a∈R Φ f0 -fa 2 emp + 1 M 2 ϕ∈S(a) log(M ) M m=1 1 0 ϕ 2 (t, X (m) )dN (m) t + log 2 (M )||ϕ|| 2 ∞   
where C is a constant depending on µ, ε, α and R.

To shed lights on this result, assume that the density of the X (m) 's is upper bounded by a constant R. In an asymptotic perspective with M → ∞, we have almost surely, M up to constants. So, we can draw the same conclusions as for the Poisson model. We have not discussed here the choice of Φ and Condition (4.3). This will be extensively done in Section 5.2 where we deal with a similar condition but in a more involved setting.

1 M M m=1 1 0 ϕ 2 (t, X (m) )dN (m) t → E 1 0 ϕ 2 (t, X (1) )f 0 (t, X (1) )Y (1) dt . But E 1 0 ϕ 2 (t, X (1) )f 0 (t, X (1) )Y (1) dt ≤ ||f 0 || ∞ E 1 0 ϕ 2 (t, X (1) )dt ≤ R||f 0 || ∞ . So, if the dictionary Φ satisfies max ϕ∈Φ ||ϕ|| ∞ = O M log M , imsart-

Applications to the case of multivariate Hawkes process

For a multivariate Hawkes model, the parameter

f 0 = (ν (m) , (h (m) ) =1,...,M ) m=1,...,M
belongs to

H = H M = f = (f (m) ) m=1,...,M | f (m) ∈ H and ||f || 2 = M m=1 ||f (m) || 2
where

H = f = (µ, (g ) =1,...,M ) | µ ∈ R , g with support in [0, 1] and ||f || 2 = µ 2 + M =1 1 0 g 2 (x)dx < ∞ .
If one defines the linear predictable transformation κ of H by

κ t (f ) = µ + M =1 t- t-1 g (t -u)dN ( ) u , (5.1)
then the transformation ψ on H is given by

ψ (m) t (f ) = κ t (f (m) ).
The first oracle inequality of Theorem 2 provides theoretical guaranties of our Lasso methodology in full generality and in particular, even if inhibition takes place (see Section 1.2.3). Since Ω V,B and Ω c are observable events, we know whether the oracle inequality holds. However we are not able to determine P(Ω V,B ) and P(Ω c ) in the general case. Therefore, in Sections 5.1 and 5.2, we assume that all interaction functions are nonnegative and that there exists f 0 in H so that for any m and any t,

λ (m) t = ψ (m) t (f 0 ).
We also assume that the process is observed on [-1, T ] with T > 1.

Some useful probabilistic results for multivariate Hawkes processes

In this paragraph, we present some particular exponential results and tail controls for Hawkes processes. As far as we know, these results are new: They constitute the generalization of [START_REF] Reynaud-Bouret | Some non asymptotic tail estimates for Hawkes processes[END_REF] to the multivariate case. In this paper, they are used to control P(Ω c c ) and P(Ω c V,B ) but they may be of independent interest. Since the functions h (m) 's are nonnegative, a cluster representation exists. We can indeed construct the Hawkes process by the Poisson cluster representation (see [START_REF] Daley | An introduction to the theory of point processes[END_REF]) as follows:

• Distribute ancestral points with marks = 1, . . . , M according to homogeneous Poisson processes with intensities ν ( ) on R. • For each ancestral point, form a cluster of descendant points. More precisely, starting with an ancestral point at time 0 of a certain type, we successively build new generations as Poisson processes with intensity h (m) (. -T ), where T is the parent of type (the corresponding children being of type m). We will be in the situation where this process becomes extinguished and we denote by H the last children of all generations, which also represents the length of the cluster. Note that the number of descendants is a multitype branching process (and there exists a branching cluster representation (see [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF][START_REF] Daley | An introduction to the theory of point processes[END_REF][START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF])) with offspring distributions being Poisson variables with means

γ ,m = 1 0 h (m) (t)dt.
The essential part we need is that the expected number of offsprings of type m from a point of type is γ ,m . With Γ = (γ ,m ) ,m=1,...,M , the theory of multitype branching processes gives that the clusters are finite almost surely if the spectral radius of Γ is strictly smaller than 1. In this case, there is a stationary version of the Hawkes process by the Poisson cluster representation. Moreover, if Γ has spectral radius strictly smaller than 1, one can provide a bound on the number of points in a cluster. We denote by P the law of the cluster whose ancestral point is of type , E is the corresponding expectation.

The following lemma is very general and holds even if the function h (m) have infinite support as long as the spectral radius Γ is strictly less than 1.

Lemma 1. If W denotes the total number of points of any type in the cluster whose ancestral point is of type then if the spectral radius of Γ is strictly smaller than 1 there exists ϑ > 0, only depending on and on Γ, such that

E (e ϑ W ) < ∞.
This easily leads to the following result, which provides the existence of the Laplace transform of the total number of points in an arbitrary bounded interval, when the functions h (m) have bounded support.

Proposition 2. Let N be a stationary multivariate Hawkes process, with compactly supported nonnegative interactions functions and such that the spectral radius of Γ is strictly smaller than 1. For any A > 0, let N [-A,0) be the total number of points of N in [-A, 0), all marks included. Then there exists a constant θ > 0, depending on the distribution of the process and on A, such that

E := E(e θN [-A,0) ) < ∞,
which implies that for all positive u

P(N [-A,0) ≥ u) ≤ Ee -θu .
Moreover one can strengthen the ergodic theorem in a non-asymptotic way, under the same assumptions. Proposition 3. Under the assumptions of Proposition 2, let A > 0 and let Z(N ) be a function depending on the points of N lying in [-A, 0). Assume that there exist b and η nonnegative constants such that

|Z(N )| ≤ b(1 + N η [-A,0)
), where N [-A,0) represents the total number of points of N in [-A, 0), all marks included. We denote S the shift operator, meaning that Z • S t (N ) depends now in the same way as Z(N ) on some points that are now the points of N lying in [t -A, t).

We assume E[|Z(N )|] < ∞ and for short, we denote

E(Z) = E[Z(N )].
Then, for any α > 0, there exists a constant T (α, η, f 0 , A) > 1 such that for T > T (α, η, f 0 , A), there exist C 1 , C 2 , C 3 and C 4 positive constants depending on α, η, A and f 0 such that

P T 0 [Z • S t (N ) -E(Z)]dt ≥ C 1 σ T log 3 (T ) + C 2 b(log(T )) 2+η ≤ C 4 T α , with σ 2 = E([Z(N ) -E(Z)] 2 1 N [-A,0) ≤ Ñ ) and Ñ = C 3 log(T ).
Finally, to deal with the control of P(Ω c ), we shall need the next result. First, we define a quadratic form Q on H by

Q(f , g) = E P (κ 1 (f )κ 1 (g)) = E P 1 T T 0 κ t (f )κ t (g)dt , f , g ∈ H. (5.2)
We have:

(m) (t) < ∞ (5.3)
then there is a constant ζ > 0 such that for any f ∈ H,

Q(f , f ) ≥ ζ||f || 2 .
We are now ready to establish oracle inequalities for multivariate Hawkes processes.

Lasso for Hawkes processes

In the sequel, we still consider the main assumptions of the previous subsection: We deal with a stationary Hawkes process whose intensity is given by (1.3) such that the spectral radius of Γ is strictly smaller than 1 and (5.3) is satisfied. We recall that the components of Γ are the γ ,m 's with

γ ,m = 1 0 h (m) (t)dt.
One of the main results of this section is to link properties of the dictionary (mainly orthonormality but also more involved assumptions) to properties of G (the control of Ω c ). To do so, let us define for all f ∈ H,

||f || ∞ = max max m=1,...,M |µ (m) |, max m, =1,...,M ||g (m) || ∞ .
Then, let us set ||Φ|| ∞ := max{||ϕ|| ∞ , ϕ ∈ Φ}. The next result is based on the probabilistic results of Section 5.1.

Proposition 5. Assume that the Hawkes process is stationary, that (5.3) is satisfied and that the spectral radius of Γ is strictly smaller than 1. Let r Φ be the spectral radius of the matrix H defined by

H = m |µ (m) ϕ ||µ (m) ρ | + M =1 1 0 |(g ϕ ) (m) ||(g ρ ) (m) |(u)du ϕ,ρ∈Φ
.

Assume that Φ is orthonormal and that

A Φ (T ) := r Φ ||Φ|| 2 ∞ |Φ|[log(||Φ|| ∞ ) + log(|Φ|)] log 5 (T ) T → 0 (5.4)
when T → ∞. Then, for any β > 0, there exists C 1 > 0 depending on β and f 0 such that with c = C 1 T , we have

P(Ω c c ) = O(|Φ| 2 T -β ).
Up to logarithmic terms, (5.4) is similar to (4.3) with M replaced with T . The dictionary Φ can be built via a dictionary (Υ k ) k=1,...,K of functions of L 2 ([0, 1]) (that may depend on T ) in the following way. A function ϕ = (µ

(m) ϕ , ((g ϕ ) (m) ) ) m belongs to Φ if
and only if only one of its M + M 2 components is non zero and in this case,

• if µ (m) ϕ = 0, then µ (m) ϕ = 1, • if (g ϕ ) (m) = 0, then there exists k ∈ {1, . . . , K} such that (g ϕ ) (m) = Υ k .
Note that |Φ| = M + KM 2 . Furthermore, assume from now on that (Υ k ) k=1,...,K is orthonormal in L 2 ([0, 1]). Then Φ is also orthonormal in H endowed with ||.||. Before going further, let us discuss Assumption (5.4). First note that the matrix H is block diagonal. The first block is the identity matrix of size M . The other M 2 blocks are identical to the matrix:

H K = |Υ k1 (u)||Υ k2 (u)|du 1≤k1,k2≤K
.

So, if we denote rK the spectral radius of H K , we have:

r Φ = max(1, rK ).
We analyze the behavior of rK with respect to K. Note that for any k 1 and any k 2 ,

(H K ) k1,k2 ≥ 0.
Therefore, rK ≤ sup

||x|| 1 =1 ||H K x|| 1 ≤ max k1 k2 (H K ) k1,k2 .
We now distinguish three types of orthonormal dictionaries (remember that M is viewed as a constant here):

• Let us consider regular histograms. The basis is composed of the functions

Υ k = δ -1/2 1 ((k-1)δ,kδ] with Kδ = 1. Therefore ||Φ|| ∞ = δ -1/2 = √ K.
But H K is the identity matrix and rK = 1. Hence (5.4) is satisfied as soon as

K 2 log(K) log 5 (T ) T → 0 when T → ∞, which is satisfied if K = o √ T log 3 (T ) . • Assume that ||Φ|| ∞ is
bounded by an absolute constant (Fourier dictionaries satisfy this assumption). Since rK ≤ K, (5.4) is satisfied as soon as

K 2 log(K) log 5 (T ) T → 0 when T → ∞, which is satisfied if K = o √ T log 3 (T ) . • Assume that (Υ k ) k=1,.
..,K is a compactly supported wavelet dictionary where resolution levels belong to the set {0, 1, . . . , J}. In this case, K is of the same order as 2 J , ||Φ|| ∞ is of the same order as 2 J/2 and it can be established that rK

≤ C2 J/2
where C is a constant only depending on the choice of the wavelet system (see [START_REF] Härdle | Wavelets, Approximation and Statistical Applications[END_REF] for further details). Then, (5.4) is satisfied as soon as To apply Theorem 2, it remains to control Ω V,B . Note that

K 5/2 log(K) log 5 (T ) T → 0 when T → ∞, which is satisfied if K = o T 2/
ψ (m) t (ϕ) = 1 if µ (m) ϕ = 1 t- t-1 Υ k (t -u)dN ( ) u if (g ϕ ) (m) = Υ k .
Let us define

Ω N = for all t ∈ [0, T ], for all m ∈ {1, . . . , M } we have N (m) [t-1,t] ≤ N .
We therefore set

B ϕ = 1 if µ (m) ϕ = 1 and B ϕ = ||Υ k || ∞ N if (g ϕ ) (m) = Υ k . (5.5)
Note that on Ω N , for any ϕ ∈ Φ,

sup t∈[0,T ],m |ψ (m) t (ϕ)| ≤ B ϕ .
Now, for each ϕ ∈ Φ, let us determine V ϕ that constitutes an upper bound of

M ϕ = M m=1 T 0 [ψ (m) t (ϕ)] 2 dN (m) t
.

Note that only one term in this sum is non-zero. We set

V ϕ = T N if µ (m) ϕ = 1 and V ϕ = ||Υ k || 2 ∞ T N 3 if (g ϕ ) (m) = Υ k , (5.6)
where T denotes the smallest integer larger than T . With this choice, one has that Ω N ⊂ Ω V,B , which leads to the following result.

Corollary 3. Assume that the Hawkes process is stationary, that (5.3) is satisfied and that the spectral radius of Γ is strictly smaller than 1. With the choices (5.5) and (5.6),

P(Ω V,B ) ≥ P(Ω N ) ≥ 1 -C 1 T exp(-C 2 N ),
where C 1 and C 2 are positive constants depending on f 0 .

If N log(T ), then for all β > 0,

P(Ω c V,B ) ≤ P(Ω c N ) = o(T -β ).
We are now ready to apply Theorem 2.

Corollary 4. Assume that the Hawkes process is stationary, that (5.3) is satisfied and that the spectral radius of Γ is strictly smaller than 1. Assume that the dictionary Φ is built as previously from an orthonormal family (Υ k ) k=1,...,K . With the notations of imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

Theorem 2, let B ϕ be defined by (5.5) and d ϕ be defined accordingly with x = α log(T ).

Then, with probability larger than

1 -4(M + M 2 K)   log 1 + µ T N α log(T ) log(1 + ε) + 1   T -α -P(Ω c N ) -P(Ω c c ), 1 T || f -f 0 || 2 T ≤ C inf a∈R Φ    1 T f 0 -f a 2 T + ϕ∈S(a) log(T )(ψ(ϕ)) 2 • N T T 2 + B 2 ϕ log 2 (T ) T 2    ,
where C is a constant depending on f 0 , µ, ε, and α.

From an asymptotic point of view, if the dictionary also satisfies (5.4), and if N = log 2 (T ) in (5.5), then for T large enough with probability larger than 1 -

C 1 K log(T )T -α 1 T || f -f 0 || 2 T ≤ C 2 inf a∈R Φ    1 T f 0 -f a 2 T + log 3 (T ) T ϕ∈S(a) 1 T ||ϕ|| 2 T + log 7/2 (T ) √ T ||Φ|| 2 ∞    ,
where C 1 and C 2 are constants depending on M , f 0 , µ, ε, and α.

We express the oracle inequality by using 1 T ||.|| T simply because, when T goes to +∞, by ergodicity of the process (see for instance [START_REF] Daley | An introduction to the theory of point processes[END_REF], and Proposition 3 for a non asymptotic statement),

1 T ||f || 2 T = M m=1 1 T T 0 (κ t (f (m) )) 2 dt -→ M m=1 Q(f (m) , f (m) )
under assumptions of Proposition 5. Note that the right hand side is a true norm on H by Proposition 4. Note also that

log 7/2 (T ) √ T ||Φ|| 2 ∞ T →∞ → 0,
as soon as (5.4) is satisfied for the Fourier basis and compactly supported wavelets. It is also the case for histograms as soon as K = o √ T log 7/2 (T ) . Therefore, this term can be viewed as a residual one. In those cases, the last inequality can be rewritten as

1 T || f -f 0 || 2 T ≤ C inf a∈R Φ    1 T f 0 -f a 2 T + log 3 (T ) T ϕ∈S(a) 1 T ||ϕ|| 2 T    ,
for a different constant C, the probability of this event tending to 1 as soon as α ≥ 1/2 in the Fourier and histogram cases and α ≥ 2/5 for the compactly supported wavelet basis. Once again, as mentioned for the Poisson or Aalen models, the right hand side corresponds to a classical "bias-variance" trade off and we obtain a classical oracle inequality up to the logarithmic terms. Note that asymptotics is now with respect to T and not with respect to M as for Poisson or Aalen models. So, the same result, namely Theorem 2, allows to consider both asymptotics.

Simulations for the multivariate Hawkes process

This section is devoted to illustrations of our procedure on simulated data of multivariate Hawkes processes and comparisons with the well-known adaptive Lasso procedure proposed by [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. We consider the general case and we do no longer assume that the functions h (m) are nonnegative as in Section 5. However, if the parameter ν (m) is large with respect to the h (m) 's, then ψ (m) (f 0 ) is nonnegative with large probability and therefore λ (m) = ψ (m) (f 0 ) with large probability. Hence, Theorem 2 implies that f is close to f 0 .

Description of the Data

As mentioned in the introduction, Hawkes processes can be used in neuroscience to model the action potentials of individual neurons. So, we perform simulations whose parameters are close, to some extent, to real neuronal data. For a given neuron m ∈ {1, . . . , M }, we recall that its activity is modeled by a point process N (m) whose intensity is

λ (m) t = ν (m) + M =1 t- -∞ h (m) (t -u)dN ( ) (u) + .
The interaction function h (m) represents the influence of the past activity of the neuron on the neuron m. The spontaneous rate ν (m) may somehow represent the external excitation linked to all the other neurons that are not recorded. It is consequently of crucial importance not only to correctly infer the interaction functions, but also to reconstruct the spontaneous rates accurately. Usually, activity up to 10 neurons can be recorded in a "stationary" phase during a few seconds (sometimes up to one minute). Typically, the points frequency is of the order of 10-80 Hz and the interaction range between points is of the order of a few milliseconds (up to 20 or 40 ms). We first lead three experiments in the pure excitation case where all the interaction functions are nonnegative by simulating multivariate Hawkes processes (two with M = 2, one with M = 8) based on these typical values. More precisely, we take for any m ∈ {1, . . . , M }, ν (m) = 20 and the interaction functions h (m) are defined as follows (supports of all the functions are assumed to lie in the interval [0, 0.04]):

• Experiment 1: M = 2 and piecewise constant functions.

h (1) 1 = 30 × 1 (0,0.02] , h (1) 2 = 30 × 1 (0,0.01] , h (2) 1 = 30 × 1 (0.01,0.02] , h (2)
2 = 0. In this case, each neuron depends on the other one. The spectral radius of the matrix Γ is 0.725. • Experiment 2: M = 2 and "smooth" functions. In this experiment, h (2)

1 (x) = 1 0.008 √ 2π e -(x-0.02) 2 2 * 0.004 2 × 1 (0,0.04] (x), h (2)
2 (x) = 0.

In this case, each neuron depends on the other one as well. The spectral radius of the matrix Γ is 0.711. • Experiment 3: M = 8 and piecewise constant functions.

h (1) 2 = h (1) 3 = h (2) 2 = h (3) 1 = h (3) 2 = h (5) 8 = h (6) 5 = h (7) 6 = h (8) 7 = 25 × 1 (0,0.02]
and all the other 55 interaction functions are equal to 0. Note in particular that this leads to 3 independent groups of dependent neurons {1, 2, 3}, {4} and {5, 6, 7, 8}. The spectral radius of the matrix Γ is 0.5.

We also lead one experiment in the pure inhibition case where all the interaction functions are nonpositive:

• Experiment 4: M = 2. In this experiment, the interaction functions are the opposite of the functions introduced in Experiment 2. We take for any m ∈ {1, . . . , M }, ν (m) = 60 so that ψ t (f 0 ) is positive with high probability.

For each simulation, we let the process "warm up" during 1 second to reach the stationary state1 . Then the data are collected by taking recordings during the next T seconds. For instance, we record about 100 points per neuron when T = 2 and 1000 points when T = 20. Figure 1 shows two instances of data sets with T = 2.

Description of the methods

To avoid approximation errors when computing the matrix G, we focus on a dictionary (Υ k ) k=1,...,K whose functions are piecewise constant. More precisely, we take Υ k = δ -1/2 1 ((k-1)δ,kδ] with δ = 0.04/K and K, the size of the dictionary, is chosen later.

Our practical procedure strongly relies on the theoretical one based on the d ϕ 's defined in (2.8), with x, µ and ε to be specified. First, using Corollary 4, we naturally take x = α log(T ). Then, three hyperparameters would need to be tuned, namely α, µ and ε, if we directly used the Lasso estimate of Theorem 2. So, for simplifications, we implement our procedure by replacing the Lasso parameters d ϕ with

dϕ (γ) = 2γ log(T )(ψ(ϕ)) 2 • N T + γ log(T ) 3 sup t∈[0,T ],m |ψ (m) t (ϕ)|,
where γ is a constant to be tuned. Besides taking x = α log(T ), our modification consists in neglecting the linear part

B 2 ϕ x µ-φ(µ) in V µ and replacing B ϕ with sup t∈[0,T ],m |ψ (m) t
(ϕ)|. Then, note that, up to these modifications, the choice γ = 1 corresponds to the limit case where α → 1, ε → 0 and µ → 0 in the definition of the d ϕ 's (see the comments after Theorem 2). Note also that, under the slight abuse consisting in identifying B ϕ with sup t∈[0,T ],m |ψ (m) t (ϕ)|, for every parameter µ, ε and α of Theorem 2 with x = α ln(T ), one can find two parameters γ and γ such that dϕ (γ) ≤ d ϕ ≤ dϕ (γ ). Therefore, this practical choice is consistent with the theory and tuning hyperparameters reduces to only tuning γ. Our simulation study will provide sound answers to the question of tuning γ.

We compute the Lasso estimate by using the shooting method of [START_REF] Fu | Penalized regressions: the bridge versus the Lasso[END_REF] and the R-package Lassoshooting. In particular, we need to invert the matrix G. In all simulations, this matrix was invertible, which is consistent with the fact that Ω c happens with large probability. Note also that the value of c, namely the smallest eigenvalue of G, can be very small (about 10 -4 ) whereas the largest eigenvalue is potentially as large as 10 5 , both values highly depending on the simulation and on T . Fortunately, those values are not needed to compute our Lasso estimate. Since it is based on Bernstein type inequalities, our Lasso method is denoted B in the sequel.

Due to their soft thresholding nature, Lasso methods are known to underestimate the coefficients [START_REF] Meinshausen | Relaxed Lasso[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. To overcome biases in estimation due to shrinkage, we propose a two steps procedure, as usually suggested in the literature: Once the support of the vector has been estimated by B, we compute the ordinary least-square estimator among the vectors a having the same support, which provides the final estimate. This method is denoted BO in the sequel.

Another popular method is adaptive Lasso proposed by Zou [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. This method overcomes the flaws of standard Lasso by taking 1 -weights of the form where p > 0, γ > 0 and âo ϕ is a preliminary consistent estimate of the true coefficient. Even if the shapes of the weights are different, the latter are data-driven and this method constitutes a natural competitive method with ours. The most usual choice, which is adopted in the sequel, consists in taking p = 1 and the ordinary least squares estimate for the preliminary estimate (see [START_REF] Huang | Adaptive Lasso for sparse highdimensional regression models[END_REF][START_REF] Van De Geer | The adaptive and the thresholded Lasso for potentially misspecified models (and a lower bound for the Lasso)[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]). Then, penalization is stronger for coefficients that are preliminary estimated by small values of the ordinary least square estimate. In the literature, the parameter γ of adaptive Lasso is usually tuned by cross-validation, but this does not make sense for Hawkes data that are fully dependent. Therefore, a preliminary study has been performed to provide meaningful values for γ. Results are given in the next section. This adaptive Lasso method is denoted A in the sequel and AO when combined with ordinary least squares in the same way as for BO.

d a ϕ (γ) = γ 2|â o ϕ | p , imsart-
Simulations are performed in R. The computational time is small (merely a few seconds for one estimate even when M = 8, T = 20 and K = 8 on a classical laptop computer), which constitutes a clear improvement with respect to existing adaptive methods for Hawkes processes. For instance, the "Islands" method2 of [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF] is limited to the estimation of one or two dozens of coefficients at most, because of an extreme computational memory cost whereas here when M = 8 and K = 8, we can easily deal with M + KM 2 = 520 coefficients.

Results

First, we provide in Figure 2 reconstructions by using the OLS estimate on the whole dictionary, which corresponds to the case where all the weights d ϕ are null. As expected, reconstructions are not sparse and also bad due to a small signal to noise ratio (remember that T = 2). Now let us consider methods leading to sparsity. A precise study over 100 simulations has been carried out corresponding to Experiments 1 and 3 for which we can precisely check if the support of the vector â is the correct one. For each method, we have selected 3 values for the hyperparameter γ based on results of preliminary simulations. Before studying mean squared errors, we investigate the following problems that are stated in order of importance. We wonder whether our procedure can identify:

the dependency groups. Recall that two neurons belong to the same group if and only if they are connected directly or through the intermediary of one or several neurons. This issue is essential from the neurobiological point of view since knowing interactions between two neurons is of capital importance. -the non-zero interaction functions h (m) 's and non-zero spontaneous rates ν (m) 's. For , m ∈ {1, . . . , M }, the neuron has a significative direct interaction on neuron m if and only if h (m) = 0;

the non-zero coefficients of non-zero interaction functions. This issue is more mathematical. However, it may provide information about the maximal range for direct interactions between two given neurons or about the favored delay of interaction.

Note that the dependency groups are the only features that can be detected by classical analysis tools of neuroscience, such as the Unitary Events method [START_REF] Grün | Detecting unitary events without discretization in time[END_REF]. In particular, to the best of our knowledge, identification of the non-zero interaction functions inside a dependency group is a problem that has not been solved yet as far as we know. Results for our method and for adaptive Lasso can be found in Table 1. This preliminary study also provides answers for tuning issues. The line "DG" gives the number of correct identifications of dependency groups. For instance, for M = 8, "DG" gives the number of simulations for which the 3 dependency groups {1, 2, 3}, {4} and {5, 6, 7, 8} are recovered by the methods. When M = 2, both methods correctly find that neurons 1 and 2 are dependent, even if T = 2. When 8 neurons are considered, the estimates given for T = 2 (left) and T = 20 (right). "DG" gives the number of correct identifications of dependency groups over 100 runs. "S" gives the median number of non-zero spontaneous rate estimates, "*" means that all the spontaneous rate estimates are non-zero over all the simulations. "F+" gives the median number of additional non-zero interaction functions w.r.t. the truth. "F-" gives the median number of missing non-zero interaction functions w.r.t. the truth. "Coeff+" and "Coeff-" are defined in the same way for the coefficients. "SpontMSE" is the Mean Square

Error for the spontaneous rates with or without the additional "ordinary least squares step". "InterMSE" is the analog for the interaction functions. In red, we give the optimal values.

imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013 should find 3 dependency groups. We see that even with T = 2, our method with γ = 1 correctly guesses the dependency groups for 32% of the simulations. It's close or equal to 100% when T = 20 with γ = 1 or γ = 2. The adaptive Lasso has to take γ = 1000 for T = 2 and T = 20 to obtain as convincing results. Clearly, smaller choices of γ for adaptive Lasso leads to bad estimations of the dependency groups. Next, let us focus on the detection of non-zero spontaneous rates. Whatever the experiment and the parameter γ, our method is optimal whereas adaptive Lasso misses some non-zero spontaneous rates when T = 2. Under this criterion, for adaptive Lasso, the choice γ = 1000 is clearly bad when T = 2 (the optimal value of S is S = 2 when M = 2 and S = 8 when M = 8) on both experiments, whereas γ = 2 or γ = 200 is better. Not surprisingly, the number of additional non-zero functions and additional non-zero coefficients decreases when T grows and when γ grows, whatever the method whereas the number of missing functions or coefficients increases. We can conclude from these facts and from further analysis of Table 1 that the choice γ = 0.5 for our method and the choice γ = 2 for the adaptive Lasso are wrong choices of the tuning parameters. In conclusion of this preliminary study, our method with γ = 1 or γ = 2 seems a good choice and is robust with respect to T . When T = 20, the optimal choice for adaptive Lasso is γ = 1000. When T = 2, the choice is not so clear and depends on the criterion we wish to favor. Now let us look at mean squared errors (MSE). Since the spontaneous rates do not behave like the other coefficients, we split the MSE in two parts: one for the spontaneous rates:

SpontMSE = M m=1 (ν (m) -ν (m) ) 2 ,
and one for interactions:

InterMSE = M m=1 M =1 ( ĥ(m) (t) -h (m) (t)) 2 dt.
We still report the results for B, BO, A and AO in Table 1. Our comments mostly focus on cases where the results for the previous study are good. First, note that results on such cases are better by using the second step (OLS). Furthermore, MSE is increasing with γ for B and A, since underestimation is stronger when γ increases. This phenomenon does not appear for two step procedures, which leads to a more stable MSE. For adaptive Lasso, when T = 2, the choice γ = 200 leads to good MSE, but the MSE are smaller for BO with γ = 1. When T = 20, the choice γ = 1000 for AO leads to results that are of the same magnitude as the ones obtained by BO with γ = 1 or 2. Still for T = 20, results for the estimate B are worse than results for A. It is due to the fact that shrinkage is larger in our method for the coefficients we want to keep than shrinkage of adaptive Lasso that becomes negligible as soon as the true coefficients are large enough. However the second step overcomes this problem. Note also that a more thorough study of the tuning parameter γ has been performed by [START_REF] Bertin | Adaptive Dantzig density estimation[END_REF] where it is mathematically proved that the choice γ < 1 leads to very degenerate estimates in the density setting. Their method for choosing Lasso parameters being analogous to ours, it seems coherent to obtain worse MSE for γ = 0.5 than for γ = 1 or γ = 2, at least for BO. The boundary γ = 1 in their simulation study seems to be a robust choice there, and it seems to be the case here too. We now provide some reconstructions by using Lasso methods. Figures 3 and4 give the reconstructions corresponding to Experiment 2 (M = 2) with K = 8 for T = 2 and T = 20 respectively. The reconstructions are quite satisfying. Of course, the quality improves when T grows. We also note improvements by using BO instead of B. For adaptive Lasso, improvements by using the second step are not significative and this is the reason why we do not represent reconstructions with AO. Graphs of the right hand side of Figure 3 illustrate the difficulties of adaptive Lasso to recover the exact support of interactions functions, namely h the sake of clarity, we only represent reconstructions for the first 4 neurons. From the estimation point of view, this illustration provides a clear hierarchy between the methods: BO seems to achieve the best results and B the worst. Finally, Figure 6 shows that even in the inhibition case, we are able to recover the negative interactions.

Conclusions

With respect to the problem of tuning our methodology based on Bernstein type inequalities, our simulation study is coherent with theoretical aspects since we achieve our best results by taking γ = 1, which constitutes the limit case of assumptions of Theorem 2.

For practical aspects, we recommend the choice γ = 1 even if γ = 2 is acceptable. More importantly, this choice is robust with respect to the duration of recordings, which is not the case for adaptive Lasso. Implemented with γ = 1, our method outperforms adaptive Lasso and it is able to recover the dependency groups, the non-zero spontaneous rates, the non-zero functions and even the non-zero coefficients as soon as T is large enough. Most of the time, the two step procedure BO seems to achieve the best results for parameter estimation.

It is important to note that the question of tuning adaptive Lasso remains open. Some values of γ allow us to obtain very good results but they are not robust with respect to T , which may constitute a serious problem for practitioners. In the standard regression setting, this problem may be overcome by using cross-validation on independent data, which somehow estimates random fluctuations. But in this multivariate Hawkes setup, independence assumptions on data cannot be made and this explains the problems for tuning adaptive Lasso. Our method based on Bernstein type concentration inequalities takes into account those fluctuations. It also takes into account the nature of the coefficients and the variability of their estimates which differ for spontaneous rates on the one hand and coefficients of interaction functions on the other hand. The shape of weights of adaptive Lasso does not incorporate this difference, which explains the contradictions for tuning the method when T = 2. For instance, in some cases, adaptive Lasso tends to estimate some spontaneous rates to zero in order to achieve better performance on the imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013 interaction functions.

Proofs

This section is devoted to the proofs of the results of the paper. Throughout, C is a constant whose value may change from line to line.

Proof of Theorem 1

The proof of Theorem 1 is standard (see for instance [START_REF] Bunea | Sparsity Oracle Inequalities for the Lasso[END_REF]), but for the sake of completeness, we give it. We use ||.|| 2 for the Euclidian norm of R Φ . Given a recall that

f a = ϕ∈Φ a ϕ. Then, we have f = f â, a b = ψ(f a ) • N T and a Ga = ||f a || 2 T = ||ψ(f a )|| 2 proc . Then, -2ψ(f â) • N T + ||f â|| 2 T + 2d |â| ≤ -2ψ(f a ) • N T + ||f a || 2 T + 2d |a|. So, ||ψ(f â) -λ|| 2 proc = ||ψ(f â)|| 2 proc + ||λ|| 2 proc -2 < ψ(f â), λ >proc ≤ ||ψ(fa)|| 2 proc + ||λ|| 2 proc + 2ψ(f â -fa) • NT +2d (|a| -|â|) -2 < ψ(f â), λ >proc = ||ψ(fa) -λ|| 2 proc + 2 < ψ(fa -f â), λ >proc +2ψ(f â -fa) • NT + 2d (|a| -|â|) = ||ψ(fa) -λ|| 2 proc + 2ψ(fa -f â) • (Λ -N )T + 2d (|a| -|â|) = ||ψ(fa) -λ|| 2 proc + 2 ϕ∈Φ (aϕ -âϕ)ψ(ϕ) • (Λ -N )T + 2d (|a| -|â|) ≤ ||ψ(fa) -λ|| 2 proc + 2 ϕ∈Φ |aϕ -âϕ| × | bϕ -bϕ| + 2d (|a| -|â|) .
Using (2.5), we obtain:

||ψ(f â) -λ|| 2 proc ≤ ||ψ(f a ) -λ|| 2 proc + 2 ϕ∈Φ d ϕ |a ϕ -âϕ | + 2 ϕ∈Φ d ϕ (|a ϕ | -|â ϕ |) ≤ ||ψ(f a ) -λ|| 2 proc + 2 ϕ∈Φ d ϕ (|a ϕ -âϕ | + |a ϕ | -|â ϕ |) .
||ψ(f â) -λ|| 2 proc ≤ ||ψ(f a ) -λ|| 2 proc + 2 ϕ∈S(a) d ϕ (|a ϕ -âϕ | + |a ϕ | -|â ϕ |) ≤ ||ψ(f a ) -λ|| 2 proc + 4 ϕ∈S(a) d ϕ (|a ϕ -âϕ |) ≤ ||ψ(f a ) -λ|| 2 proc + 4||â -a|| 2   ϕ∈S(a) d 2 ϕ   1/2
.

We now use the assumption on the Gram matrix given by (2.4) and the triangular inequality for ||.|| T , which yields

||â -a|| 2 2 ≤ c -1 (â -a) G (â -a) = c -1 ||f â -f a || 2 T ≤ 2c -1 ||ψ(f â) -λ|| 2 proc + ||ψ(f a ) -λ|| 2 proc .
Let us take α ∈ (0, 1). Since for any x ∈ R and any y ∈ R, 2xy ≤ αx 2 + α -1 y 2 , we obtain:

||ψ(f â) -λ|| 2 proc ≤ ||ψ(fa) -λ|| 2 proc +4 √ 2c -1/2 ||ψ(f â) -λ|| 2 proc + ||ψ(fa) -λ|| 2 proc   ϕ∈S(a) d 2 ϕ   1/2 ≤ ||ψ(fa) -λ|| 2 proc +α ||ψ(f â) -λ|| 2 proc + ||ψ(fa) -λ|| 2 proc + 8α -1 c -1 ϕ∈S(a) d 2 ϕ ≤ (1 -α) -1   (1 + α)||ψ(fa) -λ|| 2 proc + 8α -1 c -1 ϕ∈S(a) d 2 ϕ   .
The theorem is proved just by taking an arbitrary absolute value for α ∈ (0, 1).

Proof of Theorem 2

Let us first define

T = {t ≥ 0 / sup m |ψ (m) t (ϕ)| > B ϕ }. (7.1)
Let us define the stopping time τ = inf T and the predictible process H by

H (m) t = ψ (m) t (ϕ)1 t≤τ .
Let us apply Theorem 3 to this choice of H with τ = T and B = B ϕ . The choice of v and w will be given later on. To apply this result, we need to check that for all t and all imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

ξ ∈ (0, 3), m t 0 e ξ H (m) s Bϕ λ (m)
s ds is a.s. finite. But if t > τ , then

t 0 e ξ H (m) s Bϕ λ (m) s ds = τ 0 e ξ H (m) s Bϕ λ (m) s ds + t τ λ (m) s ds,
where the second part is obviously finite (it is just Λ

(m) t -Λ (m)
τ .) Hence it remains to prove that for all t ≤ τ , | ≤ B ϕ . Since we are integrating with respect to the Lebesgue measure, the fact that it eventually does not hold in t is not a problem and

t 0 e ξ H (m) s Bϕ λ (m) s ds ≤ e ξ Λ (m) t ,
which is obviously finite a.s. The same reasoning can be applied to show that a.s. exp(ξH 2 /B 2 ) • Λ t < ∞. We can also apply Theorem 3 to -H in the same way. We obtain at the end that for all ε > 0

P |H • (N -Λ)T | ≥ 2(1 + ε) V µ x + Bϕx 3 and w ≤ V µ ≤ v and sup m,t≤T |H (m) t | ≤ Bϕ ≤ 4 log(v/w) log(1 + ε) + 1 e -x . (7.2)
But on Ω V,B it is clear that ∀t ∈ [0, T ], t ∈ T . Therefore τ ≥ T . Therefore for all t ≤ T , one also has t ≤ τ and H

(m) t = ψ (m) t (ϕ). Consequently, on Ω V,B , H • (N -Λ) T = b ϕ -bϕ and V µ = V µ ϕ .
Moreover, on Ω V,B , one has that

B 2 ϕ x µ -φ(µ) ≤ V µ ϕ ≤ µ µ -φ(µ) V ϕ + B 2 ϕ x µ -φ(µ)
.

So, we take w and v as respectively the left and right hand side of the previous inequality. Finally note that on Ω V,B , sup

m,t≤T |H (m) t | = sup m,t≤T |ψ (m) t (ϕ)| ≤ B ϕ .
Hence, we can rewrite (7.2) as follows Apply this to all ϕ ∈ Φ, we obtain that

P |b ϕ -bϕ | ≥ 2(1 + ε) V µ ϕ x + B ϕ x 3 and Ω V,B ≤ 4   log 1 + µVϕ B 2 ϕ x log(1 + ε) + 1   e -x .
P ∃ ϕ ∈ Φ s.t. |b ϕ -bϕ | ≥ d ϕ and Ω V,B ≤ 4 ϕ∈Φ   log 1 + µVϕ B 2 ϕ x log(1 + ε) + 1   e -x .
Now on the event Ω c ∩ Ω V,B ∩ {∀ϕ ∈ Φ, |b ϕ -bϕ | ≤ d ϕ }, one can apply Theorem 1. To obtain Theorem 2, it remains to bound the probability of the complementary event by

P(Ω c c ) + P(Ω c V,B ) + P ∃ ϕ ∈ Φ s.t. |b ϕ -bϕ | ≥ d ϕ and Ω V,B .

Proof of Theorem 3

First, replacing H with H/B, we can always assume that B = 1. Next, let us fix for the moment ξ ∈ (0, 3). If one assumes that almost surely for all t > 0,

M m=1 t 0 e ξH (m) s λ (m)
s ds < ∞ (i.e. that the process e ξH • Λ is well defined) then one can apply Theorem 2 of [8, p165], stating that the process (E t ) t≥0 defined for all t by

E t = exp(ξH • (N -Λ) t -φ(ξH) • Λ t )
is a supermartingale. It is also the case for E t∧τ if τ is a bounded stopping time. Hence for any ξ ∈ (0, 3) and for any x > 0, one has that So, for every ξ ∈ (0, 3), we obtain:

P(E t∧τ > e x ) ≤ e -x E(E t∧τ ) ≤ e -x ,
P M τ ≥ ξ -1 φ(ξ)H 2 • Λ τ + ξ -1 x and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x . (7.4)
Now let us focus on the event H 2 • Λ τ ≤ v where v is a deterministic quantity. We have that consequently

P M τ ≥ ξ -1 φ(ξ)v + ξ -1 x and H 2 • Λ τ ≤ v and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x .
imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

It remains to choose ξ such that ξ -1 φ(ξ)v + ξ -1 x is minimal. But this expression has no simple form. However, since 0 < ξ < 3, one can bound φ(ξ) by ξ 2 (1 -ξ/3) -1 /2. Hence we can start with

P M τ ≥ ξ 2(1 -ξ/3) H 2 • Λ τ + ξ -1 x and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x (7.5)
and also

P M τ ≥ ξ 2(1 -ξ/3) v + ξ -1 x and H 2 • Λ τ ≤ v and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x . (7.6) It remains now to minimize ξ -→ ξ 2(1-ξ/3) v + ξ -1 x.
Lemma 2. Let a, b and x be positive constants and let us consider on (0, 1/b),

g(ξ) = aξ (1 -bξ) + x ξ .
Then min ξ∈(0,1/b) g(ξ) = 2 √ ax + bx and the minimum is achieved in ξ(a, b, x)

= xb- √ ax xb 2 -a .
Proof. The limits of g in 0 + and (1/b) -are +∞. The derivative is given by

g (ξ) = a (1 -bξ) 2 - x ξ 2
which is null in ξ(a, b, x) (remark that the other solution of the polynomial does not lie in (0, 1/b)). Finally it remains to evaluate the quantity in ξ(a, b, x) to obtain the result. Now, we apply (7.6) with ξ(v/2, 1/3, x) and we obtain this well known formula which can be found in [START_REF] Shorack | Empirical processes with applications to statistics[END_REF] for instance:

P M τ ≥ √ 2vx + x/3 and H 2 • Λ τ ≤ v and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x . (7.7)
Now we would like first to replace v by its random version H 2 • Λ τ . Let w, v be some positive constants and let us concentrate on the event

w ≤ H 2 • Λ τ ≤ v. (7.8)
For all ε > 0 we introduce K a positive integer depending on ε, v and w such that

(1 + ε) K w ≥ v. Note that K = log(v/w)/ log(1 + ε) is a possible choice. Let us denote v 0 = w, v 1 = (1 + ε)w, ..., v K = (1 + ε) K w.
For any 0 < ξ < 3 and any k in {0, ..., K -1}, one has, by applying (7.5),

P Mτ ≥ ξ 2(1 -ξ/3) H 2 • Λτ + ξ -1 x and v k ≤ H 2 • Λτ ≤ v k+1 and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x .
imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

This implies that

P Mτ ≥ ξ 2(1 -ξ/3) v k+1 + ξ -1 x and v k ≤ H 2 • Λτ ≤ v k+1 and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x .
Using Lemma 2, with ξ = ξ(v k+1 /2, 1/3, x), this gives

P M τ ≥ 2v k+1 x + x/3 and v k ≤ H 2 • Λ τ ≤ v k+1 and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x . But if v k ≤ H 2 • Λ τ , v k+1 ≤ (1 + ε)v k ≤ (1 + ε)H 2 • Λ τ , so P Mτ ≥ 2(1 + ε) (H 2 • Λτ ) x + x/3 and v k ≤ H 2 • Λτ ≤ v k+1 and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x .
Finally summing on k, this gives

P M τ ≥ 2(1 + ε)(H 2 • Λ τ )x + x/3 and w ≤ H 2 • Λ τ ≤ v and sup s≤τ,m |H (m) s | ≤ 1 ≤ Ke -x . (7.9)
This leads to the following result that has interest per se. with respect to some given filtration. Let B > 0. Let H = (H (m) ) m=1,...,M be a multivariate predictable process such that for all ξ ∈ (0, 3), e ξH/B •Λ t < ∞ a.s. for all t. Let us consider the martingale defined for all t by

M t = H • (N -Λ) t .
Let v > w be positive constants and let τ be a bounded stopping time. Then for any ε, x > 0

P Mτ ≥ 2(1 + ε)(H 2 • Λτ )x + Bx 3 and w ≤ H 2 • Λτ ≤ v and sup m,t≤τ |H (m) t | ≤ B ≤ log(v/w) log(1 + ε) + 1 e -x . (7.10)
Next, we would like to replace H 2 • Λ τ , the quadratic characteristic of M , with its estimator H 2 • N τ , i.e. the quadratic variation of M . For this purpose, let us consider W t = -H 2 • (N -Λ) t which is still a martingale since the -(H (m) s ) 2 's are still predictable processes. We apply (7.4) with µ instead of ξ, noticing that on the event

{sup s≤τ,m |H (m) s | ≤ 1}, one has that H 4 • Λ τ ≤ H 2 • Λ τ . This gives that P H 2 • Λ τ ≥ H 2 • N τ + {φ(µ)/µ}H 2 • Λ τ + x/µ and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x ,
imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: [START_REF] Chagny | Adaptive warped kernel estimators[END_REF] which means that

P H 2 • Λ τ ≥ V µ and sup s≤τ,m |H (m) s | ≤ 1 ≤ e -x .
(7.11)

So we use again (7.5) combined with (7.11) to obtain that for all ξ ∈ (0, 3)

P M τ ≥ ξ 2(1 -ξ/3) V µ + ξ -1 x and sup s≤τ,m |H (m) s | ≤ 1 ≤ P M τ ≥ ξ 2(1 -ξ/3) V µ + ξ -1 x and sup s≤τ,m |H (m) s | ≤ 1 and H 2 • Λ τ ≤ V µ + + P H 2 • Λ τ ≥ V µ and sup s≤τ,m |H (m) s | ≤ 1 ≤ 2e -x .
This new inequality replaces (7.5) and it remains to replace H 2 • Λ τ by V µ in the peeling arguments to obtain as before that

P M τ ≥ 2(1 + ε) V µ x + x/3 and w ≤ V µ ≤ v and sup s≤τ,m |H (m) s | ≤ 1 ≤ 2Ke -x .
(7.12)

7.4. Proofs of the probabilistic results for Hawkes processes Thus, φ (θ) is the log-Laplace transform of the distribution of K(1) given that there is a single initial ancestral point of type . We define the vector φ(θ) by φ(θ) = (φ 1 (θ), ..., φ M (θ)). Note that φ only depends on the law of the number of children per parent, i.e. it only depends on Γ. Then

E e θ T W (n) = E e θ T W (n-1) E e θ T K(n) | K(n -1), . . . , .K(1) = E e θ T W (n-1) e φ(θ) T K(n-1)
= E e (θ+φ(θ)) T K(n-1)+θ T W (n-2)

Defining g(θ) = θ + φ(θ) we arrive by recursion at the formula

E e θ T W (n) = E e g •(n-1) (θ) T K(1)+θ T W (0)
= e φ(g where for any n, g

•(n) = g • • • • • g n times.
Or, in other words, we have the following

representation log E e θ T W (n) = g •n (θ)
of the log-Laplace transform of W (n). Below we show that φ is a contraction in a neighborhood containing 0, that is, for some r > 0 and a constant C < 1 (and a suitable norm), ||φ(s

)|| ≤ C||s|| for ||s|| ≤ r. If θ is chosen such that ||θ|| 1 -C ≤ r
we have ||θ|| ≤ r, and if we assume that g •k (θ) ∈ B(0, r) for k = 1, . . . , n -1 then

||g •n (θ)|| ≤ ||θ|| + ||φ(g •(n-1) (θ))|| ≤ ||θ|| + C||g •(n-1) (θ)|| ≤ ||θ|| 1 + C + C 2 + . . . + C n ≤ r
Thus, by induction, g •n (θ) ∈ B(0, r) for all n ≥ 1. Since n → W m (n) is increasing and goes to W m (∞) for n → ∞, with W m (∞) the total number of points in a cluster of type m, and since

W = m W m (∞) = 1 T W (∞), we have by monotone convergence that for ϑ ∈ R log E e ϑW = lim n→∞ g •n (ϑ1) .
By the previous result, the right hand side is bounded if |ϑ| is sufficiently small. This completes the proof up to proving that φ is a contraction.

To this end we note that φ is continuously differentiable (on R M in fact, but a neighborhood around 0 suffices) with derivative Dφ(0) = Γ at 0. Since the spectral radius of Γ is strictly less than 1 there is a C < 1 and, by the Householder theorem, a norm || • || on R M such that for the induced operator norm of Γ we have Note that we have not at all used the explicit formula for φ above, which is obtainable and simple since the offspring distributions are Poisson. The only thing we needed was the fact that φ is defined in a neighborhood around 0, thus that the offspring distributions are sufficiently light-tailed.

Proof of Proposition 2

We use the cluster representation, and we note that any cluster with ancestral point in [-n -1, -n] must have at least n + 1 -A points in the cluster if any of the points are to fall in [-A, 0). This follows from the assumption that all the h

H n (ϑ ) := E e ϑ max{W -n+ A ,0} ≤ P (W ≤ n -A ) + e -ϑ (n-A ) E e ϑ W ,
which is finite for |ϑ | sufficiently small according to Lemma 1. Then we can compute an upper bound on the Laplace transform of ÑA, :

Ee ϑ ÑA, ≤ n E An k=1 E e ϑ max{W n,k -n+ A ,0} | A n ≤ n EH n (ϑ ) An = n e ν (Hn(ϑ )-1)
= e ν n (Hn(ϑ )-1)

Since H n (ϑ ) -1 ≤ e -ϑ (n-A ) E e ϑ W we have n (H n (ϑ ) -1) < ∞, which shows that the upper bound is finite. To complete the proof, observe that N [-A,0) = ÑA, where ÑA, for = 1, . . . , M are independent. Since all variables are positive, it is sufficient to take θ = min ϑ .

Proof of Proposition 3

In this paragraph, the notation simply denotes a generic positive absolute constant that may change from line to line. The notation θ1,θ2,... denotes a positive constant depending on θ 1 , θ 2 , . . . that may change from line to line.

Let u = C 1 σ log 3/2 (T ) √ T + C 2 b(log(T )) 2+η , (7.13)
where the choices of C 1 and C 2 will be given later. For any positive integer k such that x := T /(2k) > A, we have by stationarity:

P T 0 [Z • St(N ) -E(Z)]dt ≥ u = P k-1 q=0 2qx+x 2qx [Z • St(N ) -E(Z)]dt + 2qx+2x 2qx+x [Z • St(N ) -E(Z)]dt ≥ u ≤ 2P k-1 q=0 2qx+x 2qx [Z • St(N ) -E(Z)]dt ≥ u 2 .
Similarly to [START_REF] Reynaud-Bouret | Some non asymptotic tail estimates for Hawkes processes[END_REF], we introduce ( M x q ) q a sequence of independent Hawkes processes, each being stationary with intensities per mark given by ψ (m) t

. For each q, we then introduce M x q the truncated process associated with M x q , where truncation means that we only consider the points lying in [2qx -A, 2qx + x]. So, if we set

F q = 2qx+x 2qx [Z • S t (M x q ) -E(Z)]dt, P T 0 [Z • St(N ) -E(Z)]dt ≥ u ≤ 2P k-1 q=0 Fq ≥ u 2 + 2kP Te > T 2k -A , (7.14)
where T e represents the time to extinction of the process. More precisely T e is the last point of the process if in the cluster representation only ancestral points before 0 are appearing. For more details, see section 3 of [START_REF] Reynaud-Bouret | Some non asymptotic tail estimates for Hawkes processes[END_REF]. So, denoting a l the ancestral points with marks l and H l a l the length of the corresponding cluster whose origin is a l , we have:

T e = max l∈{1,...,M } max a l a l + H l a l .
But, for any a > 0,

P(T e ≤ a) = E M l=1 a l E 1 {a l +H l a l ≤a} |a l = E M l=1 a l exp log P(H l 0 ≤ a -a l ) = E M l=1 exp 0 -∞ log(P(H l 0 ≤ a -x))d Ñ (l)
x , imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

We have P(Ω c ) ≤ q P(Ω c q ). Each Ω q can also be easily controlled. Indeed it is sufficient to split [2qx -A, 2qx + x] in intervals of size A (there are about α,A,f0 log(T ) of those) and require that the number of points in each subinterval is smaller than Ñ /2. By stationarity, we obtain that

P(Ω c q ) ≤ α,A,f0 log(T )P(N [-A,0) > Ñ /2).
Using Proposition 2 with u = Ñ /2 + 1/2, we obtain: P(Ω c q ) ≤ α,A,f0 log(T ) exp(α,A,f0 Ñ ) and P(Ω c ) ≤ α,A,f0 T exp(α,A,f0 Ñ ). (7.15) Note that this control holds for any positive choice of Ñ . Hence this gives also the following Lemma that will be used later. Lemma 3. For any R > 0,

P there exists t ∈ [0, T ] | M x q | [t-A,t) > R) ≤ α,A,f0 T exp(-α,A,f0 R).
Hence by taking Ñ = C 3 log(T ) for C 3 large enough this is smaller than α,A,f0 T -α , where α = max(α, 2).

It remains to obtain the rate of D := P( q F q ≥ u/2 and Ω). For any positive constant θ that will be chosen later, we have:

D ≤ e -θu 2 E e θ q Fq q 1 Ωq ≤ e -θu 2
q E e θFq 1 Ωq (7.16) since the variables (M x q ) q are independent. But

E e θFq 1 Ωq = 1 + θE(F q 1 Ωq ) + j≥2 θ j j! E(F j q 1 Ωq ) and E(F q 1 Ωq ) = E(F q ) -E(F q 1 Ω c q ) = -E(F q 1 Ω c q ). Next note that if for any integer l, l Ñ < sup t M x q | [t-A,t) ≤ (l + 1) Ñ then |F q | ≤ xb[(l + 1) η Ñ η + 1] + xE(Z).
imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

Hence, cutting Ω c q in slices of the type {l Ñ < sup t M x q [t-A,t) ≤ (l + 1) Ñ } and using Lemma 3, we obtain by taking C 3 large enough,

|E(F q 1 Ωq )| = |E(F q 1 Ω c q )| ≤ +∞ l=1 x(b[(l + 1) η Ñ η + 1] + |E(Z)|) × P(there exists t ∈ [0, T ] | {M x q | [t-A,t) } > Ñ ) ≤ α,A,f0 +∞ l=1 x(b[(l + 1) η Ñ η + 1] + |E(Z)|) log(T )e -α,A,f 0 l Ñ ≤ α,A,f0 +∞ l=1 x(b Ñ η + |E(Z)|) log(T )2 lη e -α,A,f 0 l Ñ ≤ α,η,A,f0 log 2 (T )b Ñ η e -α,A,f 0 Ñ 1 -2 η e -α,A,f 0 Ñ ≤ z 1 := α,η,A,f0 bT -α .
Note that in the previous inequalities, we have bounded

|E(Z)| by bE[N η [-A,0) ].
In the same way, one can bound

E(F j q 1 Ωq ) ≤ E(F 2 q 1 Ωq )z j-2 b , with z b := xb[ Ñ η + 1] + xE(Z) = α,η,A,f0 b log(T ) 1+η .
One can also note that by stationarity,

E(F 2 q 1 Ωq ) ≤ xE 2qx+x 2qx [Z • θ s (M x q ) -E(Z)] 2 1 {for all t,M x q | [t-A,t) ≤ Ñ } ds ≤ xE 2qx+x 2qx [Z • θ s (M x q ) -E(Z)] 2 1 {M x q | [s-A,s) ≤ Ñ } ds ≤ x 2 E([Z(N ) -E(Z)] 2 1 N [-A,0) ≤ Ñ ) ≤ z v := α,η,A,f0 (log(T )) 2 σ 2 .
Now let us go back to (7.16). We have that

D ≤ exp   - θu 2 + k ln   1 + θz 1 + j≥2 z v z j-2 b θ j j!     ≤ exp   -θ u 2 -kz 1 + k j≥2 z v z j-2 b θ j j!   ,
using that ln(1 + u) ≤ u. It is sufficient now to recognize a step of the proof of the Bernstein inequality (weak version see [41, p25]). Since kz 1 = α,η,s bT 1-α /(log(T )), one can choose α > 1, C 1 and C 2 in the definition (7.13) of u (not depending on b) such that u/2 -kz 1 ≥ √ 2kz v z + 1 3 z b z for some z = C 4 log(T ), where C 4 is a constant. Hence

D ≤ exp   -θ( 2kz v z + 1 3 z b z) + k j≥2 z v z j-2 b θ j j!   .
One can choose accordingly θ (as for the proof of the Bernstein inequality) to obtain a bound in e -z . It remains to choose C 4 large enough and only depending on α, η, A and f 0 to guarantee that D ≤ e -z ≤ α,η,A,f0 T -α . This concludes the proof of the proposition.

Proof of Proposition 4

Let Q denote a measure such that under Q the distribution of the full point process restricted to (-∞, 0] is identical to the distribution under P and such that on (0, ∞) the process consists of independent components each being a homogeneous Poisson process with rate 1. Furthermore, the Poisson processes should be independent of the process on (-∞, 0]. From Corollary 5.1.2 in [START_REF] Jacobsen | Point process theory and applications. Marked point and piecewise deterministic processes[END_REF] the likelihood process is given by

L t = exp M t - m t 0 λ (m) u du + m t 0 log λ (m) u dN (m)
u and we have for t ≥ 0 the relation

E P κ t (f ) 2 = E Q κ t (f ) 2 L t , (7.17)
where E P and E Q denote the expectation with respect to P and Q respectively. Let, furthermore, Ñ1 = N [-1,0) denote the total number of points on [-1, 0). Proposition 4 will be an easy consequence of the following lemma.

E P (1 + ε) Ñ1 < ∞ for some ε > 0 then for any f , Q(f , f ) ≥ ζ||f || 2 (7.18)
for some constant ζ > 0.

Proof. We use Hölders inequality on κ 1 (f ) where 1 p + 1 q = 1. We choose q ≥ 1 (and thus p) below to make q -1 sufficiently small. For the left hand side we have by independence of the homogeneous Poisson processes that if f = (µ, (g ) =1,...,M ),

2 p L 1 p 1 and κ 1 (f ) 2 q L -1 p 1 to get E Q κ 1 (f ) 2 ≤ E Q κ 1 (f ) 2 L 1 1 p E Q κ 1 (f ) 2 L -q p 1 1 q = Q(f , f ) 1 p E Q κ 1 (f ) 2 L 1-q 1 1 q (7.
E Q κ 1 (f ) 2 = (E Q κ 1 (f )) 2 + V Q κ 1 (f ) = µ + 1 0 g (u)du 2 + 1 0 g (u) 2 du.
Exactly as on page 32 in [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF] there exists c > 0 such that

E Q κ 1 (f ) 2 ≥ c µ 2 + 1 0 g 2 (u)du = c ||f || 2 . (7.20)
To bound the second factor on the right hand side in (7.19) we observe, by assumption, that we have the lower bound

L 1 ≥ e M (1-b) e (d-aM )N1 e -aM Ñ1
on the likelihood process. Under Q we have that (κ 1 (f ), N 1 ) and Ñ1 are independent, and with ρ = e (q-1)(aM -d) and ρ = e (q-1)(aM ) we get that

E Q κ 1 (f ) 2 L 1-q 1 ≤ e (q-1)M (b-1) E Q ρ Ñ1 E Q κ 1 (f ) 2 ρ N1 .
Here we choose q such that ρ is sufficiently close to 1 to make sure that E Q ρ Ñ1 = E P ρ Ñ1 < ∞ (see Proposition 2). Moreover, by Cauchy-Schwarz' inequality

κ 2 1 (f ) ≤ µ 2 + 1- 0 g 2 (1 -u)dN ( ) u (1 + N 1 ). ( 7

.21)

Under Q the point processes on (0, ∞) are homogeneous Poisson processes with rate 1 and N 1 , the total number of points, is Poisson. This implies that conditionally on (N

1 , . . . , N

) = (n (1) , . . . , n (M ) ) the n (m) -points for the m'th process are uniformly distributed on [0, 1], hence Proposition 2 proves that there exists ε > 0 such that E P (1 + ε) Ñ1 < ∞. This completes the proof of Proposition 4. || 2 = T ζ||f a || 2 by definition of the norm on H. Since Φ is an orthonormal system, this implies that E(a Ga) ≥ T ζ||a|| 2 . Hence, to show that Ω c is a large event for some c > 0, it is sufficient to show that for some 0 < < ζ, with high probability, for any a ∈ R Φ , |a Ga -a E(G)a| ≤ T ||a|| 2 Since for > 0, dN

E Q κ 1 (f ) 2 L 1-q 1 ≤ µ 2 + 1 0 g 2 (u)du e (q-1)M (b-1) E Q ρ Ñ1 E Q (1 + N 1 ) 2 ρ N1 c = c ||f || 2 . ( 7
c ||f || 2 ≤ (c ) 1 q ||f || 2 q Q(f , f ) 1 p or by rearranging that Q(f , f ) ≥ ζ||f || 2 with ζ = (c ) p /(c ) p-
( )
u is stationnary, one can replace E(dN It is now sufficient to take ω small enough and then T 0 large enough to obtain (7.

| + δ -1 log M ϕ 1 ∞ ϕ 2 ∞ ] ≤ D,β (δM r Φ + δ -1 Φ 2 ∞ |Φ| log M ) a 2 2 .
We choose δ =

Φ 2 ∞ |Φ| log M M rΦ
, so that with probability larger than 1 -

2|Φ| 2 M -β , 1 M |a Ga -E[a Ga]| ≤ D,β Φ 2 ∞ r Φ |Φ| log M M a 2 2 .
We use (4.1) and (4.3) to conclude as for Proposition 5 and we obtain Proposition 1.

Proof of Corollary 3

First let us cut [-1, T ] in T + 2 intervals I's of the type [a, b) such that the first T + 1 intervals are of length 1 and the last one is of length strictly smaller than 1 (eventually it is just a singleton). Then, any interval of the type [t -1, t] for t in [0, T ] is included into the union of two such intervals. Therefore the event where all the N I 's are smaller than u = N /2 is included into Ω N . It remains to control the probability of the complementary of this event. By stationarity, all the first N I 's have the same distribution and satisfy Proposition 2. The last one can also be viewed as the truncation of a stationary point process to an interval of length smaller than 1. Therefore the exponential inequality of Proposition 2 also applies to the last interval. It remains to apply T + 2 times this exponential inequality and to use a union bound.

Proof of Corollary 4

As in the proof of Proposition 3, we use the notation . The non-asymptotic part of the result is just a pure application of Theorem 2, with the choices of B ϕ and V ϕ given by (5.5) and (5.6). The next step consists in controlling the martingale ψ(ϕ) with B = B 2 ϕ and τ = T and where τ is defined in (7.1) (see the proof of Theorem 2). The assumption to be fulfilled is checked as in the proof of Theorem 2. But as previously, on Ω V,B , H • (N -Λ) T = ψ(ϕ) 2 • (N -Λ) T and also

H 2 • Λ T = ψ(ϕ) 4 • Λ T . Moreover on Ω N ⊂ Ω V,B H 2 • Λ T = ψ(ϕ) 4 • Λ T ≤ v := T M (max m ν (m) + N max m, h (m) )B 4 ϕ .
Recall that x = α log(T ). So on Ω V,B , with probability larger than 1 -(M + KM 2 )e -x = 1 -(M + KM 2 )T -α , one has that for all ϕ ∈ Φ,

ψ(ϕ) 2 • N T ≤ ψ(ϕ) 2 • Λ T + √ 2vx + B 2 ϕ x 3 .
So that for all ϕ ∈ Φ,

ψ(ϕ) 2 • N T ≤ M,f0 N ||ϕ|| 2 T + ||Φ|| 2 ∞ N 2 T N log(T ) .
Also, since N = log 2 (T ), one can apply Corollary 3, with β = α. We finally choose c as in Proposition 5. This leads to the result.

Theorem 2 .

 2 Let N = (N (m) ) m=1,...,M be a multivariate counting process with predictable intensities λ (m) t and almost surely finite corresponding compensator Λ (m) t . Define Ω V,B = for any ϕ ∈ Φ, sup t∈[0,T ],m |ψ (m) t (ϕ)| ≤ B ϕ and (ψ(ϕ)) 2 • N T ≤ V ϕ , for positive deterministic constants B ϕ and V ϕ and

Theorem 3 .

 3 Let N = (N (m) ) m=1,...,M be a multivariate counting process with predictable intensities λ (m) t and corresponding compensator Λ (m) t

  ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013 and since the distribution of

Proposition 4 .

 4 Under the assumptions of Proposition 2, if the function parameter f 0 satisfies min m∈{1,...,M } ν (m) > 0 and max l,m∈{1,...,M } sup t∈[0,1]

  h

  h

Figure 1 .

 1 Figure 1. Raster plots of two data sets with T = 2 corresponding to Experiment 2 on the left and Experiment 3 on the right. The x-axis correspond to the time of the experiment. Each line with ordinate m corresponds to the points of the process N (m) . From bottom to top, we observe 124 and 103 points for Experiment 2 and 101, 60, 117, 38, 73, 75, 86 and 86 points for Experiment 3.

  h

Figure 2 .

 2 Figure 2. Reconstructions corresponding to Experiment 2 with the OLS estimate with T = 2 and K = 8. Each line m represents the function h (m) , for = 1, 2. The spontaneous rates associated with each line m are given above the graphs where S * denotes the true spontaneous rate and its estimator is denoted by SO. The true interactions functions are plotted in black whereas the OLS estimates are plotted in magenta.

  imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date:[START_REF] Chagny | Adaptive warped kernel estimators[END_REF]

Figure 3 .

 3 Figure 3. Reconstructions corresponding to Experiment 2 with T = 2 and K = 8. Each line m represents the function h (m) , for = 1, 2. The spontaneous rates estimation associated with each line m is given above the graphs: S * denotes the true spontaneous rate and its estimators computed by using B, BO and A respectively are denoted by SB, SBO and SA. The true interactions functions (in black) are reconstructed by using B, BO and A providing reconstructions in green, red and blue respectively. We use γ = 1 for B and BO and γ = 200 for A.

for T = 2 .Figure 4 .

 24 Figure 4. Reconstructions corresponding to Experiment 2 with T = 20 and K = 8. Same convention as in Figure 3. We use γ = 1 for B and BO and γ = 1000 for A.

Figure 5 .Figure 6 .

 56 Figure 5. Reconstructions corresponding to Experiment 3 with T = 20 and K = 8 and for the first 4 neurons. Each line m represents the function h (m), for = 1, 2, 3, 4. Same convention as in Figure3. We use γ = 1 for B and BO and γ = 1000 for A.

  imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013 Now, if ϕ ∈ S(a), |a ϕ -âϕ | + |a ϕ | -|â ϕ | = 0, and

  finite. But for all s < t, s < τ and consequently s ∈ T . Therefore |H (m) s

(7. 3 )

 3 imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date:[START_REF] Chagny | Adaptive warped kernel estimators[END_REF]

  which means thatP(ξH • (N -Λ) t∧τ -φ(ξH) • Λ t∧τ > x) ≤ e -x . Therefore P(ξH • (N -Λ) t∧τ -φ(ξH) • Λ t∧τ > x and sup s≤τ,m |H (m) s | ≤ 1) ≤ e -x . But if sup s≤τ,m |H (m) s | ≤ 1,then for any ξ > 0 and any s, φ(ξH (m) s ) ≤ (H (m) s ) 2 φ(ξ).

Proposition 6 .

 6 Let N = (N (m) ) m=1,...,M be a multivariate counting process with predictable intensities λ

7. 4 . 1 .

 41 Proof of Lemma 1Let K(n) denote the vector of the number of descendants in the n'th generation from a single ancestral point of type , define K(0) = e and let W (n) = n k=0 K(k) denote the total number of points in the first n generations. Define for θ ∈ R M φ (θ) = log E e θ T K(1) .

  Since the norm is continuous and Dφ(s) is likewise there is an r > 0 such that ||Dφ(s)|| ≤ C < 1 for ||s|| ≤ r. This, in turn, implies that φ is Lipschitz continuous in the ball B(0, r) with Lipschitz constant C, and since φ(0) = 0 we get ||φ(s)|| ≤ C||s|| for ||s|| ≤ r. This ends the proof of the lemma.

Lemma 4 .

 4 If the point process is stationary under P, if e d ≤ λ (m) t ≤ a(N 1 + Ñ1 ) + b for t ∈ [0, 1] and for constants d ∈ R and a, b > 0, and if

h

  (m) (t) < ∞ for l, m = 1, . . . , M then for t ∈ [0, 1] we have e d ≤ λ (m) t ≤ a(N 1 + Ñ1 ) + b with d = log ν (m) , a = max l sup t∈[0,1]

h

  (m) (t), b = ν (m) .

7. 5 .Q

 5 Proofs of the results of Sections 4.2 and 5.2 7.5.1. Proof of Propositions 5 and 1 We first prove Proposition 5. As in the proof of Proposition 3, we use the notation . Note that for any ϕ 1 and any ϕ 2 belonging to Φ, )κ t (ϕ 2 (m) )dt and E(G ϕ1,ϕ2 ) = T M m=1 Q(ϕ 1 (m) , ϕ 2 (m) ) by using (5.2). This implies that E(a Ga) = a E(G)a = T m 4, E(a Ga) ≥ T ζ m ||f (m) a

24 ) 2 . 2 . 1 ( 1 |

 242211 [START_REF] Comte | Adaptive estimation of the conditional intensity of marker-dependent counting processes[END_REF] implies that, with high probability, for any a ∈ R Φ ,a Ga ≥ a E(G)a -T ||a|| 2 ≥ T (ζ -)||a|| 2 ,and the choice c = T (ζ -) is convenient. So, first one has to control all the coefficients of G -E(G). For all ϕ, ρ ∈ Φ, we apply Proposition 3 to imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date:[START_REF] Chagny | Adaptive warped kernel estimators[END_REF] Note that Z only depends on points lying in [-1, 0). Therefore,|Z(N )| ≤ 2M ||ϕ|| ∞ ||ρ|| ∞ 1+ N 2[-1,0) . This leads toP 1 T G ϕ,ρ -E(G ϕ,ρ ) ≥ x ϕ,ρ ≤ β,f0 T -β with x ϕ,ρ = β,f0,M [σ ϕ,ρ log 3/2 (T )T -1/2 + ||ϕ|| ∞ ||ρ|| ∞ log 4 (T )T -1 ] probability larger than 1β,f0 |Φ| 2 T -β one has that |a Ga -a E(G)a| ≤ β,f0   ϕ,ρ∈Φ |a ϕ ||a ρ |[σ ϕ,ρ log 3/2 (T )T 1/2 + ||ϕ|| ∞ ||ρ|| ∞ log 4 (T )]   .Hence, for any positive constant δ chosen later,|a Ga -a E(G)a| ≤ β,f 0 T Now let us focus on E := ϕ,ρ∈Φ |a ϕ ||a ρ | )] 2 1 N [-1,0) ≤ Ñ ) + (E[ m ψ (m) 0 (ϕ)ψ (m) 0 (ρ)]) 2 ||ϕ|| ∞ ||ρ|| ∞ with Ñ := β,f0 log(T ). Next, ) ≤ 2M ||ϕ|| ∞ ||ρ|| ∞ (1 + N 2 [-1,0) ). Hence, if N [-1,0) ≤ Ñ = β,f0 log(T ), for T large enough, ) ≤ β,M,f0 ||ϕ|| ∞ ||ρ|| ∞ log 2 )) ≤ β,M,f0 ||ϕ|| ∞ ||ρ|| ∞ log 2 (T ).imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013Hence,E ≤ β,M,f0 log 2 (T ) |) where |f | = ((|µ (m) |, (|g (m) |) =1,...,M ) m=1,...,M . |a ϕ ||(g ϕ ) (m) |(-u)dN ( ) uIf one creates artificially a process N (0) with only one point and if we decide that (g ϕ )|a ϕ ||(g ϕ ) (m) |(-u)dN ( ) uNow we apply the Cauchy-Schwarz inequality for the measure dN ( ) , which givesN [-1,0) + 1)|(gϕ) (m) |(-u)|(gρ) (m) |(-u)dN ( ) u .imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013Now let us use the fact that for every x, y ≥ 0, η, θ > 0 that will be chosen later,xy -ηe θx ≤ y θ [log(y) -log(ηθ) -1] ,with the convention that y log(y) = 0 if y = 0. Let us apply this to x = N [-1,0) + 1 andy = |(g ϕ ) (m) |(-u)|(g ρ ) (m) |(-u). We obtain that E ≤ β,M,f 0 η log 2 (T ) M m=1 ϕ,ρ∈Φ |aϕ||aρ|E (N [-1,0) + 1)e θ(N [-1,0) +1) + β,M,f 0 θ -1 log 2 (T ) (gϕ) (m) ||(gρ) (m) |(-u) log(|(gϕ) (m) ||(gρ) (m) |(-u)) -log(ηθ) -1 dN u .

uM=1 ϕ,ρ∈Φ |a ϕ ||a ρ | 1 0| 2 . 2 . 1 T

 1221 ) by f0 du. Moreover since by Proposition 2, N [-1,0) has some exponential moments there exists θ = f0 such thatE (N [-1,0) + 1)e θ(N [-1,0) +1) = f0 .With |Φ| the size of the dictionary, this leads toE ≤ β,M,f0 η|Φ| log 2 (T )||a|| 2 2 + β,M,f0 log 2 (T ) |a ϕ ||a ρ ||µ (m) ϕ ||µ (m) ρ | log(|µ (m) ϕ ||µ (m) ρ |) -log(ηθ) -1 + (g ϕ ) (m) ||(g ρ ) (m) |(u) log(|(g ϕ ) (m) ||(g ρ ) (m) |(u)) -log(ηθ) -1 du   .Consequently, using ||Φ|| ∞ and r Φ ,E ≤ β,M,f0 η|Φ| log 2 (T )||a|| 2 2 + β,M,f0 log 2 (T )r Φ [2 log(||Φ|| ∞ ) -log(ηθ) -1]||a|| 2 We choose η = |Φ| -1 and obtain that E ≤ β,M,f0 log 2 (T )r Φ [log(||Φ|| ∞ ) + log(|Φ|)]||a|| 2 Now, let us choose δ = ω/(log 2 (T )r Φ [log(||Φ|| ∞ ) + log(|Φ|)]) where ω depends only on β, M and f 0 and will be chosen later and let us go back to (7.24): |a Ga -a E(G)a| ≤ β,M,f0 ω||a|| 2 2 + β,f0,ω r Φ [log(||Φ|| ∞ ) + log(|Φ|)] ϕ,ρ∈Φ |a ϕ ||a ρ |||ϕ|| ∞ ||ρ|| ∞ log 5 (T ) T ≤ β,M,f0 ω||a|| 2 2 + β,f0,ω ||a|| 2 2 A Φ (T ).Under assumptions of Proposition 5, for T 0 large enough and T ≥ T 0 , 1 T |a Ga -a E(G)a| ≤ β,M,f0 ω||a|| 2 2 .

1 2 ϕ 1 2 ≤ D ϕ 1 ∞ ϕ 2 ∞

 11212 [START_REF] Comte | Adaptive estimation of the conditional intensity of marker-dependent counting processes[END_REF] with < ζ and Proposition 5 is proved.Arguments for the proof of Proposition 1 are similar. So we just give a brief sketch of the proof. Now, (t,X (m) )ϕ 2 (t, X (m) )dt. Let β > 0. With probability larger than 1 -2M -β , 1 M |G ϕ1,ϕ2 -E[G ϕ1,ϕ2 ]| ≤ 2βv ϕ1,ϕ2 log M M + βb ϕ1,ϕ2 log M 3M , with b ϕ1,ϕ2 = ϕ 1 ∞ ϕ 2 ∞ , v ϕ1,ϕ2 = E (t, X (m) )ϕ 2 (t, X (m) )dt |ϕ 1 |, |ϕ 2 | ,where •, • denotes the standard L 2 -scalar product. We have just used the classical Bernstein inequality combined with (4.2). So, with probability larger than 1 -2|Φ| 2 M -β , for any vector a and any δ > 0, |a Ga -E[a Ga]| ≤ D,β ϕ1,ϕ2 |a ϕ1 ||a ϕ2 |[δM |ϕ 1 |, |ϕ 2

  Some theoretical results are established in this general setting but to go further, we shall consider in Section 5 the case where the functions h (m) are nonnegative and then λ

	(m)
	t
	is a linear function of f 0 , as in Sections 1.2.1 and 1.2.2:

Table 1 .

 1 . Numerical results of both procedures over 100 runs with K = 4. Results for Experiment 1 (top) and Experiment 3 (bottom) are

  0}where A n is the number of ancestral points in [-n -1, -n] of type and W n,k is the number of points in the respective clusters. Here the A n 's and the W n,k 's are all independent, the A n 's are Poisson distributed with mean ν and the W n,k 's are i.i.d. with the same distribution as W in Lemma 1. Moreover,

		(m) -functions have
	support in [0, 1]. With ÑA, the number of points in [-A, 0) from a cluster with ancestral
	points of type we thus have the bound
		An
	ÑA, ≤	max{W n,k -n + A ,
	n	k=1

  For the Hawkes process it follows that if ν (m) > 0 and if

	sup
	t∈[0,1]

1 . imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

  2

  • (N -Λ) T on Ω V,B .To do so, let us apply (7.7) to H such that for any m,

	H t (m)	= ψ t (m)	(ϕ) 2 1 t≤τ ,

imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date:[START_REF] Chagny | Adaptive warped kernel estimators[END_REF]

Note that since the size of the support of the interaction functions is less or equal to 0.04, the "warm up" period is

times the interaction range. imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October 3, 2013

N. R. Hansen, P. Reynaud-Bouret and V. Rivoirard

This method developed for M = 1 could easily be theoretically adapted for larger values of M , but its extreme computational cost prevents us from using it in practice.imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: October

[START_REF] Bacry | Some limit theorems for Hawkes processes and application to financial statistics[END_REF] 2013

Acknowledgements:

We are very grateful to Christine Tuleau-Malot who allowed us to use her R programs simulating Hawkes processes. The research of Patricia Reynaud-Bouret and Vincent Rivoirard is partly supported by the french Agence Nationale de la Recherche (ANR 2011 BS01 010 01 projet Calibration). The authors would like to thank the anonymous Associate Editor and Referees for helpful comments and suggestions.

where Ñ (l) denotes the process associated with the ancestral points with marks l. So,

Now, by Lemma 1, there exists some ϑ l > 0, such that c l = E (e ϑ l W ) < +∞, where W is the number of points in the cluster. But if all the interaction functions have support in [0, 1], one always have that H l 0 < W . Hence

So,

So, there exists a constant C α,f0,A depending on α, A, and f 0 such that if we take k = C α,A,f0 T / log(T ) , then

In this case x = T 2k ≈ log(T ) is larger than A for T large enough (depending on A, α, f 0 ). Now, let us focus on the first term B of (7.14), where

Let us consider some Ñ where Ñ will be fixed later and let us define the measurable events

where M x q | [t-A,t) represents the set of points of M x q lying in [t-A, t). Let us also consider Ω = ∩ 1≤q≤k Ω q . Then B ≤ P q F q ≥ u/2 and Ω + P(Ω c ).

imsart-bj ver. 2013/03/06 file: LassoHawkesFinal.tex date: [START_REF] Chagny | Adaptive warped kernel estimators[END_REF]