Agathe Guilloux

Iacopo Mastromatteo

Keywords: Convex Optimization, Stochastic Gradient Descent, Monte Carlo Markov Chain, Survival Analysis, Conditional Random Fields CHAPTER IV Constrained optimization approach Hawkes processes, convex relaxation, ADMM, compressed sensing Hawkes Process, Non-parametric estimation, GMM method, Order books, Market Microstructure

Je tiens en premier lieu à exprimer ma plus profonde gratitude envers mes directeurs de thèse Stéphane Gaï as et Emmanuel Bacry. Le premier pour son enthousiasme et son large spectre de connaissances en statistique, le second pour son Øair et son recul dans l'étude des signaux temporels. Leur soutien et leur conAEance durant ces trois ans m'ont permis de mener à bien ce travail. Merci pour tout ce que vous m'avez fait découvrir et pour tout ce que j'ai appris sous votre direction ! Je remercie Manuel Gomez Rodriguez et Niels Richard Hansen pour l'intérêt qu'ils ont porté à mon travail en acceptant de rapporter ma thèse, et m'excuse à nouveau de leur avoir volé une partie de leur été. Je suis également très honoré de la présence de Nicolas Vayatis et de Vincent Rivoirard dans mon jury de thèse.

Je suis très reconnaissant envers mes co-auteurs

Résumé

Le but de cette thèse est de montrer comment certaines méthodes d'optimisation récentes permettent de résoudre des problèmes d'estimation di ciles posés par l'étude d'évènements aléatoires dans le temps. Alors que le cadre classique de l'apprentissage supervisé traite les observations comme une collection de couples indépendants de covariables et de labels, les modèles d'évènements s'intéressent aux temps d'arrivée, à valeurs continues, de ces évènements et cherchent à extraire de l'information sur la source de donnée. Ces évènements datés sont liés par la chronologie, et ne peuvent dès lors être considérés comme indépendants. Ce simple constat justiAEe l'usage d'un outil mathématique particulier, appelé processus ponctuel, pour apprendre une structure à partir de ces évènements.

Deux exemples de processus ponctuels sont étudiés dans cette thèse. Le premier est le processus ponctuel sous-jacent au modèle de Cox à risques proportionnels : son intensité conditionnelle permet de déAEnir le ratio de risque, une quantité fondamentale dans la littérature de l'analyse de survie. Le modèle de régression de Cox relie la durée avant l'apparition d'un évènement, appelé défaillance, aux covariables d'un individu. Ce modèle peut être reformulé à l'aide du cadre des processus ponctuels. Le second est le processus de Hawkes qui modélise l'impact des évènements passés sur la probabilité d'apparition d'évènements futurs. Le cas multivarié permet d'encoder une notion de causalité entre les di érentes dimensions considérées.

Cette thèse est divisée en trois parties. La première s'intéresse à un nouvel algorithme d'optimisation que nous avons développé. Il permet d'estimer le vecteur de paramètre de la régression de Cox lorsque le nombre d'observations est très important. Notre algorithme est basé sur l'algorithme SVRG et utilise une méthode MCMC pour approcher la direction de descente. Nous avons prouvé des vitesses de convergence pour notre algorithme et avons montré sa performance numérique sur des jeux de données simulées et issues du monde réel. La deuxième partie montre que la causalité au sens de Hawkes peut être estimée de manière non-paramétrique grâce aux cumulants intégrés du processus ponctuel multivarié. Nous avons développé deux méthodes d'estimation des intégrales des noyaux du processus de Hawkes, sans faire d'hypothèse sur la forme de ces noyaux. Nos méthodes sont plus rapides et plus robustes, vis-à-vis de la forme des noyaux, par rapport à l'état de l'art. Nous avons démontré la consistance statistique de la première méthode, et avons montré que la deuxième peut être réduite à un problème d'optimisation convexe. La dernière partie met en lumière les dynamiques de carnet d'ordre grâce à la première méthode d'estimation non-paramétrique introduite dans la partie précédente. Nous avons utilisé des données du marché à terme EUREX, déAEni de nouveaux modèles de carnet d'ordre (basés sur les précédents travaux de Bacry et al.) et appliqué la méthode d'estimation sur ces processus ponctuels. Les résultats obtenus sont très satisfaisants et cohérents avec une analyse économétrique. Un tel travail prouve que la méthode que nous avons développée permet d'extraire une structure à partir de données aussi complexes que celles issues de la AEnance haute-fréquence. v

Abstract

The aim of this thesis is to show how recent optimization methods help solving tough estimation problems based on the event models. While the classical framework of supervised learning treats the observations as a collection of covariate and label independent pairs, event models only focus on the arrival dates of these events and then seek to extract information about the data source. These timestamped events are ordered chronologically and can not therefore be considered independent. This simple fact justiAEes the use of a particular mathematical tool called point process to learn some structure from these events. Two examples of point processes are studied in this thesis. The AErst is the underlying point process in the Cox model with proportional hazards: its conditional intensity allows to deAEne the risk ratio, a fundamental quantity in the literature of the survival analysis. The Cox regression model links the duration before the occurrence of an event, called failure, to an individual's covariates. This model can be reformulated using the framework of point processes. The second is the Hawkes process, which models the impact of past events on the probability of future events. The multivariate case makes it possible to encode a notion of causality between the di erent dimensions considered. This thesis is divided into three parts. The AErst focuses on a new optimization algorithm we have developed. It allows to estimate the parameter vector of the Cox regression when the number of observations is very important. Our algorithm is based on the Stochastic Variance Reduced Gradient (SVRG) algorithm and uses a Monte Carlo Markov Chain (MCMC) method to approximate the descent direction. We have proved convergence rates for our algorithm and have shown its numerical performance on simulated and real world data sets. The second part shows that the Hawkes causality can be estimated in a non-parametric way by the integrated cumulants of the multivariate point process. We have developed two methods for estimating the integrals of the kernels of the Hawkes process, without making any hypothesis about the shape of these kernels. Our methods are faster and more robust, with respect to the shape of the kernel, compared to the state-of-the-art. We have demonstrated the statistical consistency of the AErst method, and have shown that the second method can be reduced to a convex optimization problem. The last part highlights the dynamics of the order book thanks to the AErst non-parametric estimation method introduced in the previous section. We used EUREX futures data, deAEned new order book models (based on previous work by Bacry et al.) and applied the estimation method on these point processes. The results obtained are very satisfactory and consistent with an econometric analysis. This work proves that the method that we have developed makes it possible to extract a structure from data as complex as those resulting from high-frequency AEnance.

vi

Introduction

The guiding principle of this thesis is to show how the arsenal of recent optimization methods can help solving challenging new estimation problems on events models. While the classical framework of supervised learning [START_REF] Hastie | Overview of supervised learning[END_REF] treat the observations as a collection of independent couples of features and labels, events models focus on arrival timestamps to extract information from the source of data. These timestamped events are chronologically ordered and can't be regarded as independent. This mere statement motivates the use of a particular mathematical object called point process [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] to learn some patterns from events. Let us begin by presenting and motivating the questions on which we want to shed some light in this thesis.

Motivations

The amount of data being digitally collected and stored is vast and expanding rapidly. The use of predictive analytics that extract value of this data, often referred as the data revolution, has been successfully applied in astronomy [START_REF] Feigelson | Big data in astronomy[END_REF], retail sales [MB + 12] and search engines [START_REF] Chen | Business intelligence and analytics: From big data to big impact[END_REF], among others. Healthcare institutions are now also relying on data to build customized and personalized treatment models using tools from survival analysis [START_REF] Murdoch | The inevitable application of big data to health care[END_REF]. Medical research often aims at uncovering relationships between the patient's covariates and the duration until a failure event (death or other adverse e ects) happen. The information that some patients did not die during the study is obviously relevant, but can't be casted in a regression problem where one would need to observe the lifetime for all patients. This has been circumvented in [START_REF] David | Regression models and life tables (with discussion)[END_REF], one of the most cited scientiAEc paper of all time [START_REF] Van Noorden | The top 100 papers[END_REF], with its proportional hazards model that is regarded as a regression that can also extract information from censored data, i.e. patients whose failure time is not observed. An estimation procedure of the parameter vector of the regression without any assumption on the baseline hazard, regarded sometimes as a nuisance parameter, was introduced in [START_REF] Cox | Partial likelihood[END_REF] and is done via the maximization of the partial likelihood of the model. Such procedure can e ciently handle high-dimensional covariates, which happens with biostatistics data, by adding penalization terms to the criterion to minimize [START_REF] Goeman | L1 penalized estimation in the cox proportional hazards model[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. However, algorithms to maximize Cox partial likelihood does not scale well when the number of patients is high, on the contrary to most algorithms that enabled the data revolution. We might thus ask ourselves the following question:

Introduction

Question 1. How to adapt Cox proportional hazards model regression parameter estimation algorithm to the large-scale setting ?

Few years before the twentieth century, the French sociologist Durkheim already argued that human societies were like biological systems in that they were made up of interrelated components [START_REF] Durkheim | Le suicide: étude de sociologie[END_REF]. Now that our technology enabled us to be remotely connected, plenty of AEelds involve networks, like social networks, information systems, marketing, epidemiology, national security, and others. A better understanding of those large real-world networks and processes that take place over them would have paramount applications in the mentioned domains [START_REF] Rodriguez | Structure and Dynamics of Di usion Networks[END_REF]. The observation of networks often reduces to noting when nodes of the network send a message, buy a product or get infected by a virus. We often observe where and when but not how and why messages are sent over a social network. Event data from multiple providers can however help uncovering the joint dynamics and revealing the underlying structure of a system. One way to recover the inØuence structure between di erent sources is to use a kind of point process named Hawkes process [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF][START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF], whose arrival rate of events depend on the past events. Hawkes processes have been succesfully applied to model the mutual inØuence between earthquakes with di erent times and magnitudes [START_REF] Ogata | Statistical models for earthquake occurrences and residual analysis for point processes[END_REF]. Namely, it encodes how an earthquake increases the occurence's probability of new earthquakes in the form of aftershocks, via the use of Hawkes kernels. Hawkes processes also enable measuring what we call Hawkes causality i.e. the average number of events of type i that are trigerred by events of type j . Hawkes process have been succesfully applied in a broad range of domains, the two main applications model interactions within social networks [BBH12, [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF][START_REF] Iwata | Discovering latent inØuence in online social activities via shared cascade poisson processes[END_REF] and AEnancial transactions [START_REF] Bacry | Hawkes processes in AEnance[END_REF]. However, usual estimation of Hawkes causality is done by making strong assumptions on the shape of the Hawkes kernels to simplify the inference algorithm [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF]. A common assumption is the monotonic decreasing shape of the kernels (exponential or power-law), meaning that an event impact is always instantly maximal, which is non-realistic since in practice there may exist a delay before the maximal impact. This leads to the following question:

Question 2. Can we retrieve Hawkes causality without parametrizing the kernel functions ?

To answer positively to the second question, we developed two new nonparametric estimation methods for Hawkes causality, faster and which scales better with a large number of nodes. In this part, we only focus on the AErst one, for which we have proved a consistency result. Since Bowsher's pioneering work [START_REF] Bowsher | Modelling security market events in continuous time: Intensity based, multivariate point process models[END_REF], who recognized the Øexibility and the simplicity of using Hawkes processes in order to model the joint dynamics of trades and mid-price changes of the NYSE, Hawkes processes have steadily gained in popularity in the domain of high frequency AEnance, see [START_REF] Bacry | Hawkes processes in AEnance[END_REF] for a review. Indeed, taking into account the irregular occurences of transaction data requires to consider it as a point process. Besides, in the AEnancial area, plenty of features that summarize empirical AEndings are already known. For instance, the Øow of trades is known to be autocorrelated and cross-correlated with price moves. Such features called stylized facts, from the economist Nicholas Kaldor [START_REF] Kaldor | A model of economic growth[END_REF] who referred to statistical trends that need to be taken into account despite a possible lack of microscopic understanding. These stylized facts can advantageously be captured using the Outline notion of Hawkes causality. Understanding the order book dynamics is one of the core question in AEnancial statistics, and previous nonparametric representations of order books with multivariate Hawkes processes were low-dimensional because of their estimation method's complexity. The nonparametric estimation of Hawkes causality introduced in the second part of this thesis is fast and robust to kernel functions' shape, and it is natural to wonder what kind of stylized facts it can uncover from order book timestamped data.

Question 3. Can we draw a more precise picture of order book Øows dynamics using Hawkes causality's nonparamatric estimation introduced in the second part ?

Outline

Each question presented above corresponds to a part of the thesis.

In Part I, we answer Question 1 by introducing a new stochastic gradient descent algorithm applied to the maximization of regularized Cox partial-likelihood, see details below. Indeed, the regularized Cox partial log-likelihood writes as a sum of subfunctions which depend on varying length sequences of observation, on the contrary to the usual empirical risk minimization framework where subfunctions depend on one observation. Classical stochastic gradient descent algorithms are less e ective in our case. We adapt the algorithm SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] [XZ14] by adding another sampling step: each subfunction's gradient is estimated using a Monte Carlo Markov Chain (MCMC). Our algorithm achieves linear convergence once the number of MCMC iterations is bigger than an explicit lower bound. We illustrate the outperformance of our algorithm on survival datasets.

Answers to Question 2 lie in Part II where we study two nonparametric estimation procedures for Hawkes causality. Both methods are based on the computation of the integrated cumulants of the Hawkes process and taking advantage of relations between the integrated cumulants and the Hawkes causality matrix. The AErst approach relies on matching the second and third order empirical integrated cumulants with their theoretical counterparts. This is done via the minimization of the squared norm of the di erence between the two terms, which can be viewed as a Generalized Method of Moments [START_REF] Hall | Generalized method of moments[END_REF]. However, the optimization problem to solve is non-convex providing thus an approximate solution to the exact initial problem. This second approach is based on the completion of the Hawkes causality matrix using the AErst and second integrated cumulants. The relaxation of the exact problem writes as a convex optimization problem which enables us to provide the exact solution of this approximate problem.

Introduction [START_REF] Reynaud-Bouret | Goodness-of-AEt tests and nonparametric adaptive estimation for spike train analysis[END_REF] focus on the estimation of the kernel functions, and prevent order book model's dimension from being too large and/or the dataset from being too heavy. Our nonparametric method only estimates kernels' integral, involves a lighter computation and then scales better with a large number of nodes or large number of events. We also show that the Hawkes causality matrix provides a very rich summary of the system interactions. It can thus be a valuable tool in understanding the underlying structure of a system with many type of events.

Let us now rapidly review the main results of this thesis.

Part I: Large-scale Cox model

Many supervised machine learning problems can be cast into the minimization of an expected loss over a data distribution. Following the empirical risk minimization principle, the expected loss is approximated by an average of losses over training data, and a major success has been achieved by exploiting the sum-structure to design e cient stochastic algorithms [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]. Such stochastic algorithms enable a very e cient extraction of value from massive data. Applying this to large-scale survival data, from biostatistics or economics, is of course of great importance.

In Chapter I, we review the recent advances in convex optimization with Stochastic Gradient Descent (SGD) algorithms, from the pioneering work of [START_REF] Robbins | A stochastic approximation method[END_REF] to the recent variants with variance reduction [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF] [XZ14] [SSZ13] [START_REF] Roux | A stochastic gradient method with an exponential convergence _rate for AEnite training sets[END_REF]. We then introduce the notion of Point Process [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] which provides key tools for modeling events i.e. timestamps and/or locations data. We AEnally introduce the Cox proportional hazards model [START_REF] David | Regression models and life tables (with discussion)[END_REF] that relates the time that passes before some event occurs to one or more covariates via the notion of hazard rate. In Chapter II, we introduce our new optimization algorithm to help AEtting large-scale Cox model.

Background on SGD algorithms, Point Processes and Cox proportional hazards model

In this chapter, we review the classic results behind Stochastic Gradient Descent algorithms and its variance reduced adaptations. We then introduce Cox proportional hazards model.

Stochastic Gradient Descent algorithms

SGD algorithms from a general distribution A variety of statistical and machine learning optimization problems writes

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = E ª [`(µ, ª)],
where f is a goodness of AEt measure depending implicitly on some observed data, h is a regularization term that imposes structure to the solution and ª is a random variable. Typically, f is a di erentiable function with a Lipschitz gradient, whereas h might be non-smooth - First-order optimization algorithms are all variants of Gradient Descent (GD), which can be traced back to Cauchy [START_REF] Cauchy | Méthode générale pour la résolution des systemes d'équations simultanées[END_REF]. Starting at some initial point µ 0 , this algorithm minimizes a di erentiable function f by iterating the following equation

µ t +1 = µ t °¥t r f (µ t ).
(

) 1
where r f (µ) stands for the gradient of f evaluated at µ and (¥ t ) is a sequence of step sizes. Stochastic Gradient Descent (SGD) algorithms focus on the case where r f is intractable or at least time-consuming to compute. Noticing that r f (µ) writes as an expectation, one idea is to approximate the gradient in the update step (1) with a Monte Carlo Markov Chain [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF]. For instance, replacing the exact gradient r f (µ) with its MCMC estimate has enabled a sig-niAEcant step forward in training Undirected Graphical Models [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF] and Restricted Boltzmann Machines [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. This AErst form of Stochastic Gradient Descent is called Contrastive Divergence in the mentionned contexts.

SGD Algorithms from the uniform distribution Most machine learning optimization problems involve a data AEtting loss function f averaged over sample points because of the empirical risk minimization principle [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. Namely, the objective function writes

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = 1 n n X i =1 f i (µ),
where n is the number of observations, and f i is the loss associated to the i th observation. In that case, instead of running MCMC to approximate r f , one uniformly samples a random integer i between 1 and n and replace r f (µ) with r f i (µ) in the update step (1). In the large-scale setting, computing r f (µ) at each update step represents the bottleneck of the minimization algorithm, and SGD helps decreasing the computation time. Assuming that the computation of each r f i (µ) costs 1, the computation of the full gradient r f (µ) costs n, meaning that SGD's update step is n times faster than GD's one.

The comparison of the convergence rates is however di erent. Consider f twice di erentiable on R d , µ-strongly-convex, meaning that eigenvalues of the Hessian matrix r 2 f (µ) are greater than µ > 0 for any µ 2 R d , and L-smooth, meaning that the same eigenvalues are smaller than L > 0. Convergence rates with other assumptions on the function f can be found in [B + 15]. We denote µ § its minimizer and deAEne the condition number as ∑ = L/µ. The convergence rate is deAEned for iterative methods as a tight upper bound of a pre-deAEned error, and is regarded as the speed at which the algorithm converges. Denoting µ t the iterate after t steps of an iterative algorithm and considering the di erence E f (µ t ) °f (µ § ) as error, Gradient Descent's convergence rate is O(e °t /∑ ), while Stochastic Gradient Descent's one is O(∑/t ). A convergence rate of the form O(e °AEt ) with AE > 0 is called linear convergence rate since the error decrease after one iteration is at worst linear. Equivalently, convergence rates can be phrased as the total complexity to reach a AExed accuracy i.e. the number of iterations after Recently, di erent works improved Stochastic Gradient Descent using variance reduction techniques from Monte Carlo methods. The idea is to add a control variate term to the descent direction to improve the bias-variance tradeo in the approximation of the real gradient r f (µ). Those variants also enjoy linear convergence rates, and then smaller complexities (

Point processes

Point process is a useful mathematical tool to describe phenomena occuring at random locations and/or times. A point process is a random element whose values are point patterns on a set S. We present here the useful results when the set S is the interval [0, T ), and points are timestamps of events; this special case is sometimes called temporal point process. The book [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] is regarded as the main reference on point processes' theory.

Every realization of a point process ª can be written as ª = P n i =1 ± t i where ± is the Dirac measure, n is an integer-valued random variable and t i 's are random elements of [0, T ). It can be equivalently represented by a counting process

N t = R t 0 ª(s)d s = P n i =1 1 {t i ∑t } .
The usual characterization of temporal point process is done via the conditional intensity function, which is deAEned as the inAEnitesimal rate at which events are expected to occur after t , given the history of N s prior to t :

∏(t |F t ) = lim h!0 P(N t +h °Nt = 1|F t ) h
, where F t is the AEltration of the process that encodes information available up to (but not including) the time t . The most simple temporal point process is the Poisson process which assumes that the events arrive at a constant rate, which corresponds to a constant instensity function ∏ t = ∏ > 0. Note that temporal point processes can also be characterized by the distribution of interevent times i.e. the duration between two consecutive events. We remind that the distribution of interevent times of a Poisson process with intensity ∏ is an exponential distribution of parameter ∏. See the Page 41 of [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] for four equivalent ways of deAEning a temporal point process.

Two examples of temporal point process are treated in this thesis. The AErst is the point process behind Cox proportional hazards model: its conditional intensity function allows to deAEne the hazard ratio, a fundamental quantity in survival analysis literature, see [START_REF] Andersen | Statistical models based on counting processes[END_REF]. The Cox regression model relates the duration before an event called failure to some covariates. This model can be reformulated in the framework of point processes [START_REF] Andersen | Statistical models based on counting processes[END_REF]. The second 1. Part I: Large-scale Cox model is the Hawkes process which models how past events increase the probability of future events. Its multivariate version enables encoding a notion of causality between the di erent nodes. We introduce below the Cox proportional hazards model, and the Hawkes processes in Part II.

Cox proportional hazards model

Survival analysis focuses on time-to-event data, such as the death in biological organisms and failure in mechanical systems, and is now widespread in a variety of domains like biometrics, econometrics and insurance. The variable we study is the waiting time until a well-deAEned event occurs, and the main goal of survival analysis is to link the covariates, or features, of a patient to its survival time T .

Following the theory of point processes, we deAEne the intensity as the conditioned probability that a patient dies immediately after t , given that he was alive before t :

∏(t ) = lim h!0 P(t ∑ T ∑ t + h|t ∑ T ) h .
The most popular approach, for some reasons explained below, is Cox proportional hazards model [START_REF] David | Regression models and life tables (with discussion)[END_REF]. The Cox model assumes a semi-parametric form for the hazard ratio at time t for the patient i , whose features are encoded in the vector x i 2 R d :

∏ i (t ) = ∏ 0 (t ) exp(x > i µ),
where ∏ 0 (t ) is a baseline hazard ratio, which can be regarded as the hazard ratio of a patient whose covariates are x = 0. One estimation approach considers ∏ 0 as a nuisance and only estimates µ via maximizing a partial likelihood [START_REF] David | Regression models and life tables (with discussion)[END_REF]. This way of estimating suits clinical studies where physicians are only interested in the e ects of the covariates encoded in x on the hazard ratio. This can be done with computing the ratio of hazard ratios from two di erent patients:

∏ i (t ) ∏ j (t ) = exp((x i °x j ) > µ)
For that reason, Cox model is said to be a proportional hazards model.

However, maximizing this partial likelihood is a hard problem when we deal with large-scale (meaning a large number of observations n) and high-dimensional (meaning large d ) data.

To tackle the high-dimensionality, sparse penalized approaches have been considered in the literature

[Tib96] [T + 97] [Goe10]
. The problem is now to minimize the negative of the partial log-likelihood f (µ) = °`(µ) with a penalization h(µ) that makes the predictor µ to become sparse and then select variables. We will discuss this approach and the di erent models in Chapter II. On the contrary, approaches to tackle the large-scale side of the problem do not yet exist.

Introduction

SVRG beyond Empirical Risk Minimization

Survival data (y

i , x i , ± i ) n pat i =1 contains, for each individual i = 1, . . . , n pat , a features vector x i 2 R d , an observed time y i 2 R + , which is a failure time if ± i = 1 or a right-censoring time if ± i = 0. If D = {i : ± i = 1}
is the set of patients for which a failure time is observed, if n = |D| is the total number of failure times, and if R i = { j : y j ∏ y i } is the index of individuals still at risk at time y i , the negative Cox partial log-likelihood writes

°`(µ) = 1 n X i 2D h °x> i µ + log ≥ X j 2R i exp(x > j µ) ¥i (2)
for parameters µ 2 R d . Each gradient of the negative log-likelihood then writes as two nested expectations: one from an uniform distribution over D, the other over a Gibbs distribution, see Chapter II for details.

Our minimization algorithm is doubly stochastic in the sense that gradient steps are done using stochastic gradient descent (SGD) with variance reduction, and the inner expectations are approximated by a Monte Carlo Markov Chain (MCMC) algorithm. We derive conditions on the MCMC number of iterations guaranteeing convergence, and obtain a linear rate of convergence under strong convexity and a sublinear rate without this assumption.

Part II: Uncover Hawkes causality without parametrization

In Chapters III and IV, we study two methods to uncover causal relationships from a multivariate point process. We focus on one approach per chapter.

Hawkes processes

In order to model the joint dynamics of several point processes (for example timestamps of messages sent by di erent users of a social network), we will consider the multidimensional Hawkes model, introduced in 1971 in [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF] and [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF], with cross-inØuences between the di erent processes. By deAEnition a family of d point processes is a multidimensional Hawkes process if the intensities of all of its components write as linear regressions over the past of the d processes:

∏ i t = µ i + D X k=1 Z t 0 ¡ i j (t °s)d N j s .
Another way to construct Hawkes processes is to consider the following population representation, see [START_REF] Hawkes | A cluster process representation of a selfexciting process[END_REF]: individuals of type i , 1 ∑ i ∑ d , arrive as a Poisson process of intensity µ i . Every individual can have children of all types and the law of the children of type i of an individual of type j who was born or migrated in t is an inhomogeneous Poisson process of intensity ¡ i j (• °t ).

Part II: Uncover Hawkes causality without parametrization

This construction is nice because it yields a natural way to deAEne and measure the causality between events in the Hawkes model, where the integrals

g i j = Z +1 0 ¡ i j (u) du ∏ 0 for 1 ∑ i , j ∑ d .
weight the directed relationships between individuals. Namely, introducing the counting function

N i √ j t
that counts the number of events of i whose direct ancestor is an event of j , we know from [START_REF] Bacry | Hawkes processes in AEnance[END_REF] that

E[d N i √ j t ] = g i j E[d N j t ] = g i j § j d t,
(3) where we introduced § i as the intensity expectation, satisfying

E[d N i t ] = § i d t.
However in practice, the Hawkes kernels are not directly measurable from the data and these measures of causality between the di erent kinds of events are thus inaccessible.

In the literature, there are main two classes of estimation procedures for Hawkes kernels: the parametric one and the nonparametric one. The AErst one assumes a parametrization of the Hawkes kernels, the most usual assumes the kernels are decaying exponential, and estimate the parameter via the maximization of the Hawkes log-likelihood, see for example [START_REF] Bacry | A generalization error bound for sparse and low-rank multivariate hawkes processes[END_REF] or [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF]. The second one is based either on the numerical resolution of Wiener-Hopf equations which links the Hawkes kernels to its correlation structure [START_REF] Bacry | Second order statistics characterization of hawkes processes and non-parametric estimation[END_REF] (or equivalently on the approximation of the Hawkes process as an Autoregressive model and the resolution of Yule-Walker equations [START_REF] Eichler | Graphical modeling for multivariate hawkes processes with nonparametric link functions[END_REF]), or on a method of moments via the minimization of the contrast function deAEned in [START_REF] Reynaud-Bouret | Goodness-of-AEt tests and nonparametric adaptive estimation for spike train analysis[END_REF].

In Chapters III and IV, we propose two new nonparametric estimation methods to infer the integrals of the kernels using only the integrated moments of the multivariate Hawkes process. For all estimation procedures mentionned above, including ours, we need the following stability condition so that the process admits a version with a stationary intensity:

Assumption 1. The spectral norm of G = [g i j ] satisAEes ||G|| < 1.

Generalized Method of Moments approach

A recent work [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] proved that the integrated cumulants of Hawkes processes can be expressed as functions of G = [g i j ], and provided the constructive method to obtain these expressions. The AErst approach we developed in this part is a moment matching method that AEts the second-order and the third-order integrated cumulants of the process. To that end, we have designed consistent estimators of the integrated AErst, second and third cumulants of the Hawkes process. Their theoretical counterparts are polynomials of R = (I °G) °1, as shown in Introduction [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF]:

§ i = d X m=1 R i m µ m C i j = d X m=1 § m R i m R j m K i jk = d X m=1 (R i m R j m C km + R i m C j m R km +C i m R j m R km °2 § m R i m R j m R km ).
Once we observe the process N t for t 2 [0, T ], we compute the empirical integrated cumulants on windows [°H T , H T ], and minimize the squared di erence L T between the theoretical cumulants and the empirical ones. We have proven the consistency of our estimator in the limit T ! 1, once the sequence (H T ) satisAEes some conditions. Our problem can be seen as a Generalized Method of Moments [START_REF] Hall | Generalized method of moments[END_REF].

To prove the consistency of the empirical integrated cumulants, we need the following assumption:

Assumption 2. The sequence of integration domain's half-length satisAEes H

T ! 1 and H 2

T /T ! 0.
We prove in Chaper III the following theorem of consistency.

∂ °1 P °°°°! T !1

G

The numerical part, on both simulated and real-world datasets, gives very satisfying results. We AErst simulated event data, using the thinning algorithm of [START_REF] Ogata | On lewis' simulation method for point processes[END_REF], with very di erent kernel shape -exponential, power law and rectangular -and recover the true value of G for each kind of kernel. Our method is, to the best of our knowledge, the most robust with respect to the shape of the kernels. We then ran our method on the 100 most cited websites of the MemeTracker database, and on AEnancial order book data: we outperformed state-of-the-art methods on MemeTracker and extracted nice and interpretable features from the AEnancial data. Let also mention that our method is signiAEcantly faster (roughly 50 times faster) since previous methods aim at estimating functions while we only focus on their integrals.

The simplicity of the method, that maps a list of list of timestamps to a causality map between the nodes, and its statistical consistency, incited us to design new point process models of order book and capture its dynamics. The features extracted using our method have very insightful economic interpretation. This is the main purpose of the Part III.

Constrained optimization approach

The previous approach based on the Generalized Method of Moments need the AErst three cumulants to obtain enough information from the data to recover the d 2 entries of G. Assuming that the matrix G has a certain structure, we can get rid of the third order cumulant and design another estimation method using only the AErst two integrated cumulants. Plus, the resulting optimization problem is convex, on the contrary to the minimization of L T above, which enables the convergence to the global minimum. The matrix we want to estimate minimize a simple criterion f convex, typically a norm, while being consistent with the AErst two empirical integrated cumulants.

We formulate our problem as the following constrained optimization problem: On the contrary to the optimization problem of the previous chapter, the problem just stated is convex. We test this procedure on numerical simulations of various Hawkes kernels and real order book data, and we show how the criterion f impact the matrices we retrieve.

min G f (G) s.t. C = (I °G) °1L(I °G> ) °1 ||G|| < 1 g i j ∏ 0 where f (G)

Introduction

A single asset 12-dimensional Hawkes order book model

As a AErst application of the procedure described in Chapter III, we consider the following 12dimensional point process, a natural extension of the 8-dimensional point process introduced in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]:

N t = (T + t , T °t , L + t , L °t ,C + t ,C °t , T a t , T b t , L a t , L b t ,C a t ,C b t )
where each dimension counts the number of events before t :

• T + (T °): upwards (downward) mid-price move triggered by a market order.

• L + (L °): upwards (downward) mid-price move triggered by a limit order.

• C + (C °): upwards (downward) mid-price move triggered by a cancel order.

• T a (T b ): market order at the ask (bid) that does not move the price.

• L a (L b ): limit order at the ask (bid) that does not move the price.

• C a (C b ): cancel order at the ask (bid) that does not move the price.

We then use the causal interpretation of Hawkes processes to interpret our solution as a measure of the causality between events. This application of the method to this new model revealed the di erent interactions that lead to the high-frequency price mean reversion, and those between liquidity takers and liquidity makers.

For instance, one observes the e ects of T + events on other events on Figure A.1 (in the AErst columnn on the left). The most relevant interactions are the T + ! L + and T + ! L °: the latter is more intense and related to the mean-reversion of the price. Indeedn when a market order consumes the liquidity available at the best ask, two main scenarios can occur for the mid-price to change again, either the consumed liquidity is replaced, reverting back the price (mean-reverting scenario, highly probable) or the price moves up again and a new best bid is created.

A multi-asset 16-dimensional Hawkes order book model

The nonparametric estimation method introduced in Chapter III allows a fast estimation for a nonparametric methodology. We then scale up the model so as to account for events on two assets simultaneously and unveil a precise structure of the high-frequency cross-asset dynamics. We consider a 16-dimensional model, made of two 8-dimensional models of the form

N t = (P + t , P °t , T a t , T b t , L a t , L b t ,C a t ,C b t )
where the dimension P + (P °) counts upwards (downward) mid-price move triggered by any order.

We compared two couples of assets that share exposure to the same risk factors. The main empirical result of this study concerned the couple (DAX, EURO STOXX) for which price

Background on SGD algorithms, Point Processes and Cox proportional hazards model 1 SGD algorithms

Objectives that are decomposable as a sum of a number of terms come up often in applied mathematics and scientiAEc computing. They are particularly prevalent in machine learning applications, where one wants to minimize the average loss function over all observations. In the last two decades research on optimisation problems with a summation structure has focused more on the stochastic approximation setting, where the summation is assumed to be over an inAEnite set of terms [NJLS09, DS09, BCN16, Bot98]. The AEnite sum case has seen a resurgence in recent years after the discovery that there exist fast stochastic incremental gradient methods whose convergence rates are better deterministic AErst order methods. We provide a survey of fast stochastic gradient methods in the later parts of this section.

DeAEnitions

In this work, we particularly focus on problems that have convex objectives. This is a major restriction, and one at the core of much of modern optimization theory. The primary reasons for targeting convex problems are their widespread use in applications and their relative ease of solving them. For convex problems, we can almost always establish theoretical results giving a practical bound on the amount of computation time required to solve a given convex problem [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]. Convex optimisation is still of interest when addressing non-convex problems though: many algorithms that were developed for convex problems, motivated by their provably fast convergence have later been applied to non-convex problems with good empirical results [START_REF] Goodfellow | Deep Learning[END_REF].

We denote r f the gradient of f , r 2 f its Hessian matrix and ||•|| the Eucliean norm. Let now deAEne some useful notions.

I. Background on SGD algorithms, Point Processes and Cox proportional hazards model

DeAEnition 1. A function f is L-smooth with L > 0 if f is di entiable and its gradient is Lipschitz continuous, that is 8µ, µ 0 2 R d , ||r f (µ) °rf (µ 0 )|| ∑ L||µ °µ0 ||.
If the function f is twice di erentiable, the deAEnition can be equivalently written:

8µ 2 R d , |eigenvalues[r 2 f (µ)]| ∑ L.
The other assumption we will sometimes make is that of strong convexity.

DeAEnition 2. A function f is µ-strongly convex if: 8µ, µ 0 2 R d , 8t 2 [0, 1], f (t µ + (1 °t )µ 0 ) ∑ t f (µ) + (1 °t ) f (µ 0 ) °t (1 °t ) µ 2 ||µ °µ0 || 2 .
If f is di erentiable, the deAEnition can be equivalently written:

8µ, µ 0 2 R d , f (µ 0 ) ∏ f (µ) + rf (µ) > (µ 0 °µ) + µ 2 ||µ 0 °µ|| 2 .
If the function f is twice di erentiable, the deAEnition can be equivalently written:

8µ 2 R d , |eigenvalues[r 2 f (µ)]| ∏ µ.
Gradient descent based algorithms can be easily extended to non-di erentiable objectives F if they write F (µ) = f (µ) + h(µ) with f convex and di erentiable, and h convex and nondi erentiable whose proximal operator is easy to compute.

DeAEnition 3. Given a convex function h, we deAEne its proximal operator as

prox h (x) = argmin y ∑ h(y) + 1 2 ||x °y|| 2 ∏ ,
which is well-deAEned because of the strict convexity of the `2-norm.

The proximal operator can be seen as a generalization of the projection. Indeed, if h = 0 on C and h = 1 on C , prox h is exactly the projection over C . The computation of the proximal operator is also an optimization problem, but when the function h is simple enough, the proximal operator has a closed form solution. Using these proximal operators, most algorithms enjoy the same theoretical convergence rates as if the objective was di erentiable (i.e. F (µ) = f (µ)).

SGD algorithms from a general distribution

A variety of statistical and machine learning optimization problems writes

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = E ª [`(µ, ª)],
where f is a goodness of AEt measure depending implicitly on some observed data, h is a regularization term that imposes structure to the solution and ª is a random variable. Typically, f is a di erentiable function with a Lipschitz gradient, whereas h might be non-smooth (typical examples include sparsity inducing penalty).

First-order optimization algorithms are all variants of Gradient Descent (GD), which can be traced back to Cauchy [START_REF] Cauchy | Méthode générale pour la résolution des systemes d'équations simultanées[END_REF]. Starting at some initial point µ 0 , this algorithm minimizes a di erentiable function by iterating steps proportional to the negative of the gradient, as explained in Algorithm 1.

Algorithm 1 Gradient Descent (GD)

initialize µ while not converged do µ √ µ °¥r f (µ)

end while return µ

Stochastic Gradient Descent (SGD) algorithms focus on the case where r f is intractable or at least time-consuming to compute. Noticing that r f (µ) writes as an expectation like f , one idea is to approximate the gradient in the update step in Algorithm 1 with a Monte Carlo Markov Chain [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF]. Replacing the exact gradient r f (µ) with its MCMC estimate is a general approach that enabled a signiAEcant step forward in training Undirected Graphical Models [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF] and Restricted Boltzmann Machines [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. This form of Stochastic Gradient Descent is called Contrastive Divergence in the mentionned context.

Approximating the gradient of an expectation, sometimes named the score function [START_REF] Cox | Theoretical statistics[END_REF], is a recurrent task for many other problems. Among them, we can cite posterior computation in variational inference [START_REF] Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF], value function and policy learning in reinforcement learning [START_REF] Peters | Policy gradient methods[END_REF], derivative pricing [START_REF] Broadie | Estimating security price derivatives using simulation[END_REF], inventory control in operation research [START_REF] Fu | Gradient estimation[END_REF] and optimal transport theory [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF].

SGD algorithms from a uniform distribution

Most machine learning optimization problems involve a data AEtting loss function f averaged over the uniform distribution, for instance when f is the average loss function over each observation of the data set. Namely, the optimization problem to solve writes

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = 1 n n X i =1 f i (µ),
where n is the number of observations, and f i is the loss associated to the i th observation. In that case, instead of running MCMC to approximate r f , one uniformly samples a random integer i between 1 and n and replace r f (µ) with r f i (µ) in the update step, as shown in Algorithm 2. In the litterature, Stochastic Gradient Descent implicitly refers to the uniform distribution case. In the large-scale setting, computing r f (µ) at each update step represents the bottleneck of the minimization algorithm, and SGD helps decreasing the computation time.

Algorithm 2 Stochastic Gradient Descent (SGD) initialize µ as the zero vector while not converged do

pick i ª U [n] µ √ µ °¥r f i (µ)

end while return µ

Assuming the computation of each r f i (µ) costs 1, the computation of the full gradient r f (µ) costs n, meaning SGD's update step is n times faster than GD's one.

The comparison of the convergence rates is however di erent. Consider f L-smooth and convex and denote µ § its minimizer. We deAEne the condition number ∑ = L/µ. The convergence rate is measured via the di erence f (µ t ) °f (µ §

). Using the algorithm Gradient Descent with ¥ = 1/L, the convergence rates are:

f (µ t ) °f (µ § ) ∑ O µ 1 t ∂ , f (µ t ) °f (µ § ) ∑ O °e°t/∑ ¢ if f is µ-strongly convex.
The latter convergence rate which geometrically decrease the error is called linear convergence rate since the error decrease after one iteration is at worst linear. The convergence (in expectation) of the sequence (µ t ) produced by the algorithm Stochastic Gradient Descent need the step sizes to decrease to zero a speciAEc way, see [START_REF] Robbins | A stochastic approximation method[END_REF] for a general characterization. The convergence rate of stochastic algorithms is measured via the di erence E f (µ t ) °f (µ §

).

Assuming each function f i is L-Lipschitz (and not L-smooth) and f is convex, denoting

µ t = 1 t P t u=1
µ u , the convergence rates of Stochastic Gradient Descent are:

E f (µ t ) °f (µ § ) ∑ O µ 1 p t ∂ with ¥ t = 1 L p t , E f (µ t ) °f (µ § ) ∑ O ≥ ∑ t ¥ with ¥ t = 1 µt if f is µ-strongly convex.
Convergence rates with other assumptions on the function f can be found in [B + 15]. Recently, di erent works improved Stochastic Gradient Descent using variance reduction techniques from Monte Carlo methods. The idea is to add a control variate term to the descent direction to improve the bias-variance tradeo in the approximation of the real gradient r f (µ). Those variants also enjoy linear convergence rates with constant step-sizes.

1. SGD algorithms

SGD with Variance Reduction

The control variable is a variance reduction technique used in Monte Carlo methods [START_REF] Glasserman | Monte Carlo methods in AEnancial engineering[END_REF]. Its principle consists in estimating the population mean E(X ) while reducing the variance of sample of X by using a sample from another variable Y with known expectation. We deAEne a family of estimators

Z AE = AE(X °Y ) + E(Y ) AE 2 [0, 1],
whose expectation and variance equal

E(Z a ) = AEE(X ) + (1 °AE)E(Y ), V(Z a ) = AE 2 [V(X ) + V(Y ) °2cov(X , Y )].
The case AE = 1 provides an unbiased estimator, while 0 < AE < 1 implies Z AE to be biased with reduced variance. This control variates is particularly useful when Y is positively correlated with X .

The authors of [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] observed that the variance induced by SGD's descent direction can only decrease to zero if decreasing step sizes are used, which prevents from linear convergence rate. In their work, they propose a variance reduction approach on the descent direction so as to use constant step sizes and obtain a linear convergence rate. The algorithms SAG [START_REF] Roux | A stochastic gradient method with an exponential convergence _rate for AEnite training sets[END_REF][START_REF] Schmidt | Minimizing AEnite sums with the stochastic average gradient[END_REF], SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF][START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF] and SDCA [START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss minimization[END_REF] can be phrased with the variance reduction approach described above. Update steps of SAG, SAGA and SVRG with i ª U [n] respectively write this way:

(SAG) µ √ µ °¥ √ r f i (µ) °yi n + 1 n n X j =1 y j ! , (SAGA) µ √ µ °¥ √ r f i (µ) °yi + 1 n n X j =1 y j ! , (SVRG) µ √ µ °¥ √ r f i (µ) °rf i ( μ) + 1 n n X j =1 r f j ( μ) ! .
From the control variate interpreation, we observe that SAG's descent direction is a biased estimate (AE = 1/n) of the gradient r f (µ), while SAGA's and SVRG's ones are unbiased (AE = 1).

Stochastic Average Gradient (SAG)

At each iteration, the algorithm SAG [START_REF] Roux | A stochastic gradient method with an exponential convergence _rate for AEnite training sets[END_REF] computes one gradient r f i with the up-to-date value of µ, like SGD, and then descend in the direction of the average of the most recently computed gradients r f j with equals weights, see Algorithm 3. Even though some gradients in the summation haven't been updated recently, the algorithm enjoys a linear convergence rate in the strongly-convex case. SAG can be regarded as a stochastic version of Incremental Average Gradient [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF], which has the same update with a di erent constant factor, and with cyclic computation of the gradient instead of randomised. The convergence rates in the convex and strongly-convex cases with ¥ = 1/(16L) respectively involves the average iterate µ t and the iterate µ t :

E f (µ t ) °f (µ § ) ∑ O µ 1 t ∂ E f (µ t ) °f (µ § ) ∑ O ≥ e °t °1 8n ^1 16∑ ¢¥ if f is µ-strongly convex.
The algorithm SAG is adaptative to the level of convexity of the problem, as it may be used with the same step size on both convex and strongly convex problems.

Algorithm 3 Stochastic Average Gradient (SAG)

initialize µ as the zero vector,

y i = rf i (µ) for each i while not converged do µ √ µ °¥ n P n j =1 y j pick i ª U [n] y i √ rf i (µ)
end while return µ

Stochastic Variance Reduced Gradient (SVRG)

The SVRG algorithm [XZ14, JZ13] is a recent stochastic gradient algorithm with variance reduction with linear convergence rate, given in Algorithm 4. Unlike SAG and SAGA, there is another parameter m to tune, which controls the update frequency of the control variate μ. The algorithm S2GD [START_REF] Konecn | Semi-stochastic gradient descent methods[END_REF] was developed at the same time, and has the same update as SVRG. The di erence lies in the update of the control variate μ:

• Option I: μ is the average of the µ values from the last m iterations, used in [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF].

• Option II: μ is a randomly sampled µ from the last m iterations, used for S2GD [START_REF] Konecn | Semi-stochastic gradient descent methods[END_REF].

Consider f µ-strongly convex, a step size ¥ < 1/(2L), and assume m is su cently large so that

Ω = 1 µ¥(1 °2L¥)m + 2L¥ 1 °2L¥ < 1,
then the SVRG algorithm has a linear convergence rate if t is a multiple of m:

E f ( μt ) °f (µ § ) ∑ O °Ωt/m ¢ .
Let us mention that SVRG does not require the storage of full gradients, on the contrary to SDCA, SAG and SAGA. The algorithm just stores the gradient r f ( μ) and re-evaluates the gradient r f i ( μ) at each iteration.

SGD algorithms

Algorithm 4 Stochastic Variance Reduced Gradient (SVRG) initialize µ and μ as zero vectors, t as zero while not converged do

pick i ª U [n] µ √ µ °¥(r f i (µ) °rf i ( μ) + rf ( μ)) t √ t + 1 if t is a multiple of m then
update μ with option I or II end if end while return µ SAGA The algorithm SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], described in Algorithm 5, enjoys a linear convergence rate in the strongly convex case, like SAG and SVRG, but it has the advantage with respect to SAG that it allows non-smooth penalty terms such as `1 regularization. The proof of the convergence rate is easier as well, especially because SAG's descent direction is a biased estimate of the gradient, while SAGA's one is unbiased. As SAG, the algorithm SAGA maintains the current iterate µ and a table of historical gradients.

The convergence rate of the algorithm SAGA writes:

E f (µ t ) °f (µ § ) ∑ O ≥ n t ¥ with ¥ = 1 3L , E||µ t °µ § || 2 ∑ O ≥ e °t 2(n+∑) ¥ with ¥ = 1 2(µn + L) if f is µ-strongly convex.
Algorithm 5 SAGA initialize µ as the zero vector,

y i = rf i (µ) for each i while not converged do pick i ª U [n] µ √ µ °¥ ≥ r f i (µ) °yi + 1 n P n j =1 y j ¥ y i √ rf i (µ)

end while return µ

Composite case In the paragraphs above, we gave the convergence rates of the algorithm in the smooth case i.e. when the objective function to minimize is a smooth function. When the objective function is not smooth, one writes it as the sum of its smooth part f (µ) and its non-smooth part h(µ). One can easily adapt the previous algorithms by computing the gradient of the smooth part f and then project the iterate using the proximal operator of the non-smooth part h. This adds a projection step µ √ prox h (µ) at the end of each iteration.

Point Processes

Point processes are useful to describe phenomena occurring at random locations and/or times.

A point process is a random element whose values are point patterns on a set S. We present here the deAEnitions and the useful results from point processes' theory. For further details, the book [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] is regarded as the main reference in the area of point processes.

DeAEnitions

Let S be a locally compact metric space equipped with its Borel ae-algebra B. Let X S be the set of locally AEnite counting measures on S, and N S the smallest ae-algebra on X S such that all point counts f B : X S ! N, ! 7 ! #(! \ B ) are measurable for B relatively compact in B. A point process on S is a measurable map ª from a probability space (≠, F , P) to the measurable space (X S , N S ).

Every realization of a point process ª can be written as ª =

P n i =1 ± X i
where ± is the Dirac measure, n is an integer-valued random variable and X i 's are random elements of S. A point process can be equivalenty represented by a counting process:

N (B ) := R B ª(x)d x,

Temporal Point Processes

A particular interesting case of point processes is given when S is the time interval [0, T ), which we will call a temporal point process. Here, a realization is simply a set of time points:

ª = P n i =1 ± t i .
With a slight notation abuse we will write ª = {t 1 ,..., t n } where each t i is a random time before T , and we deAEne N t = P ø2ª 1 ø∑t the associated counting process. The conditional intensity function is the usual way to characterize temporal point processes where the present depends on the past. It is deAEned as the expected inAEnitesimal rate at which events are expected to occur after t , given the history of the counting process N t prior to t .

Point Processes

Namely,

∏(t |F t ) = lim d t!0 P(N t +d t °Nt = 1|F t ) d t
, where F t is the natural AEltration of the process, it represents the information available up to (but not including) the time t . The conditional intensity function is sometimes denoted ∏ § (t ). The most simple temporal point process is the homogeneous Poisson process which assumes that the events arrive at a constant rate, which corresponds to a constant intensity function

∏(t |F t ) = ∏ § (t ) = ∏ > 0.
More generally, we deAEne the inhomogeneous Poisson process for which the conditional intensity function depends on t but not on the history i.e.

∏(t |F t ) = ∏ § (t ) = ∏(t ).
The conditional intensity turns out to be interesting for multiple reasons. First, it is a convenient characterization of a temporal point process since it describes what is locally happening at t and is easy to interpret as an instantaneous probability. Secondly, the conditional intensity can be used for simulating a temporal point process: the basic idea is to simulate a Poisson process and use the cumulative conditional intensity to time scale the interevent times [START_REF] Ogata | On lewis' simulation method for point processes[END_REF]. Thirdly, the likelihood function can be expressed on closed form using the conditional intensity: if the point process is deAEned on [0, T ), then the likelihood and the log-likelihood functions are given by

L(ª) = √ n Y i =1 ∏ § (t i ) ! exp µ °ZT 0 ∏ § (s)d s ∂ , log L(ª) = n X i =1 log ∏ § (t i ) °ZT 0 ∏ § (s)d s.
Finally, the conditional intensity function is useful for many other purposes, like a goodnessof-AEt test known as residual analysis for point processes [START_REF] Ogata | Statistical models for earthquake occurrences and residual analysis for point processes[END_REF], or the conditional distribution of interevent times between events [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF]. We can also deAEne the compensator §(t ) of the point process, with respect to F t , as the integral of the conditional intensity function:

§(t ) = R t 0 ∏ § (s)d s. We remind that N t ° §(t ) is then a F t -martingale.
We remind that the distribution of interevent times of a Poisson process with intensity ∏ is an exponential distribution of parameter ∏. More generally, we denote f § (t ) the conditional probability density function of the interevent time, t n the last event that occured and T the random next one, F § (t ) = P(t n ∑ T ∑ t |F t ) the conditional cumulative density function, and

S § (t ) = 1 °F § (t ) = P(T ∏ t |F t ) the survival function. Now, ∏ § (t ) = lim h!0 P(t ∑ T ∑ t + h|T ∏ t ) h = lim h!0 1 h P(t ∑ T ∑ t + h) P(T ∏ t ) = lim h!0 µ 1 h f § (t )h S § (t ) + o(1) ∂ = f § (t ) 1 °F § (t )
.

I. Background on SGD algorithms, Point Processes and Cox proportional hazards model

Conversely, we can write the likelihood function of the next event using the conditional intensity function:

f § (t ) = ∏ § (t ) exp µ °Zt t n ∏ § (s)d s ∂ .
This last formula enables writing a point process's realization's likelihood, already introduced above.

3 Cox proportional hazards model

Survival analysis

Survival analysis focuses on time-to-event data, such as the death in biological organisms and failure in mechanical systems, and is now widespread in a variety of domains like biometrics, econometrics and insurance [START_REF] Andersen | Statistical models based on counting processes[END_REF]. The variable we study is the waiting time until a well-deAEned event occurs, and the main goal of survival analysis is to link the covariates, or features, of a patient to its survival time. We denote T the random variable of the time of death, we deAEne the survival function as:

S(t ) = P(t ∑ T ).
However, and fortunately, not all a ected patients die during a medical study and some patients can also leave the study before its end: we say that these observations are right-censored, in the sense that for some units the event of interest has not occured at the time the data are analyzed. The information about censored individual is incomplete, but it is still an information because one knows that an individual survived at least until the date he left the study.

We will only study this kind of censoring in this part1 . Let us now consider the probabilistic formulation for our framework: let T be a non-negative random variable representing the waiting time until the occurrence of an event (we will refer to this event as failure and to this waiting time as failure time). However, we don't always observe the random variable T since the patient can leave the study -before its death -at time C called the censoring time.

Actually, we do observe T ^C and we know if the patient died or left the study i.e. we know

± =
{T ∑C } . We also assume that T and C are independent. We can now describe the model using counting processes.

Let (≠, F , P) be a probability space and (F t ) t ∏0 a AEltration satisfying the usual conditions. Let N be a point process with compensator § with respect to (F

t ) t ∏0 so that N ° § is a (F t )
t ∏0 -martingale. We denote (T 1 , ..., T n ) i .i .d . copies of the random variable of interest T , corresponding to n di erent patients of a medical study for instance, (C 1 , ...,C n ) i .i .d . copies of the censoring variable C and we deAEne for each patient i :

± i = {T i ∑C i } , the counting process N i (t ) = ± i {T i ^Ci ∑t } , and Y i (t ) = {T i ^Ci ∏t } ,
which is a predictable process. To understand the behavior of the counting process N i (t ), we introduce its intensity AE i (t ) deAEned as 3. Cox proportional hazards model the conditional probability that the patient i dies immediately after t , given that he was alive before t :

AE(t ) = lim h!0 P(t ∑ T ∑ t + h|t ∑ T ) h = °S0 (t ) S(t )
Since the process can jump only once, the intensity of N i (t ) takes the form AE

i (t ) = ∏ i (t )Y i (t ),
where ∏ i (t ) is called the hazard ratio. We also introduce the cumulative hazard § i (t ) = R t 0 ∏ i (s)d s, which can be seen as the sum of the risks faced from 0 to t . Survival analysis generally aims at estimating either S(t ) or ∏(t ) (or §(t )) given the observations of n individuals. Many approachs exist: the parametric one, which assumes that the functions can be described with a AEnite and small number of parameters, the nonparametric one, which assumes that the function of interest belongs to a certain class of smooth functions and the semi-parametric one, that has parametric and non-parametric components. The most popular approach, for some reasons explained below, is Cox proportional hazards model. The Cox model [START_REF] David | Regression models and life tables (with discussion)[END_REF] assumes a semi-parametric form for the hazard ratio at time t for the patient i , whose features are encoded in the vector x i 2 R d :

∏ i (t ) = ∏ 0 (t ) exp(x > i µ)
where ∏ 0 (t ) is a baseline hazard ratio, which can be regarded as the hazard ratio of a patient whose covariates x = 0. Two estimation approachs exist: either estimating ∏ 0 and µ which can be done via maximizing the full likelihood of the model [START_REF] Ren | Full likelihood inferences in the cox model: an empirical likelihood approach[END_REF] [She15], or considering ∏ 0 a nuisance and only estimating µ via maximizing a partial likelihood L(µ) [START_REF] David | Regression models and life tables (with discussion)[END_REF]. This way of estimating suits clinical studies where physicians are only interested in the e ects of the covariates encoded in x on the hazard ratio. This can be done with computing the ratio of hazard ratios from two di erent patients:

∏ i (t ) ∏ j (t ) = exp((x i °x j ) > µ)
For that reason, Cox model is said to be a proportional hazards model.

However, maximizing such functions is a hard problem when we deal with large-scale (meaning large n) and high-dimensional (meaning large d ) data. To tackle to high-dimensionality, sparse penalized approaches have been considered in the literature [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] [

T + 97] [Goe10].
The problem is now to minimize the negative of the partial log-likelihood °`(µ) with a penalization that make the predictor µ become sparse and then select variables. We will discuss further this approach and the di erent models. On the contrary, approaches to tackle the large-scale side of the problem do not yet exist. We give an answer to this question in the following chapter.

Existing methods

The maximization of the partial likelihood L P (µ) introduced in [START_REF] David | Regression models and life tables (with discussion)[END_REF] enables the estimation of µ -without the estimation of ∏ 0 . The partial likelihood writes:

L P (µ) = n Y i =1 √ exp(x > i µ) P j 2R i exp(x > j µ) ! ± i (1)

I. Background on SGD algorithms, Point Processes and Cox proportional hazards model

We prove in appendix that the negative of the partial log-likehood is convex, then the issue of AEnding the µ that match our data can be expressed as a classical convex optimization problem. We will consider the problem of maximizing the partial likelihood of the Cox model in the rest of this chapter.

However, in case of large-scale (meaning large n) and high-dimensional (meaning large p) data, this function becomes hard to maximize. To tackle to high-dimensionality, sparse penalized approaches have been considered in the literature. The problem is now to minimize the negative of the partial log-likelihood °`(µ) + pen(µ) i.e.

1 n n X i =1 ± i " °x> i µ + log √ X j 2R i exp(x > j µ) !# + pen(µ)
where pen(µ) is a penalization term that make the predictor µ become sparse and then select variables. For instance, the sparse penalties Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] [

T + 97], Elastic-Net [SFHT11] [YZ12],
SCAD [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], Adaptative Lasso [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF], Graphical Lasso [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], SLOPE [BvdBS + 15] and others.

Indeed, the Lasso penalty [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] pen lasso (µ) = ∏||µ|| 1 can be used to obtain a penalized partial likelihood estimator b µ [START_REF] Goeman | L1 penalized estimation in the cox proportional hazards model[END_REF]. The lasso penalty tends to select only a few nonzero coe cients and does not handle well very correlated predictors: it will pick one and ignore the other.

Another well-known penalty called Ridge penalty pen ridge

(µ) = ∏ 2 ||µ|| 2
2 tends to shrink all coe cients to zero and give equal weights to very correlated predictors. Zhou and Hastie [START_REF] Simon | Regularization paths for cox's proportional hazards model via coordinate descent[END_REF] combined the strengths of the two approaches with the Elastic-Net penalty, where AE 2 [0, 1] controls the behavior of the penalty: pen e-net

(µ) = ∏ µ AE||µ|| 1 + 1 2 (1 °AE)||µ|| 2 2 ∂
The authors of [START_REF] Gopakumar | Stabilizing sparse cox model using clinical structures in electronic medical records[END_REF] studied electronic medical records and used a sparse penalty which encodes the a priori relationship between predictors i and j : A i j = 1 if predictors i and j share a temporal or well-known relation, A i j = 0 otherwise.

pen(µ) = ∏ 1 ||µ|| 1 + 1 2 ∏ 2 X i , j A i j (µ i °µj ) 2
These methods handle the high-dimensional side of the dataset, but don't look relevant when the number of patients n is large. Indeed, the higher the number of examples n, the higher the time to compute the sum of loss functions (here, the negative of the penalized log-likelihood) and time can be the limiting factor when one envisions very large datasets.

In the next chapter, we introduce a new stochastic algorithm with variance reduction that enables a faster minimization of the negative partial likelihood of the Cox model.

CHAPTER II

Large-scale Cox model

Abstract

We introduce a doubly stochastic proximal gradient algorithm for optimizing a AEnite average of smooth convex functions, whose gradients depend on numerically expensive expectations. Indeed, the e ectiveness of SGD-like algorithms relies on the assumption that the computation of a subfunction's gradient is cheap compared to the computation of the total function's gradient. This is true in the Empirical Risk Minimization (ERM) setting, but can be false when each subfunction depends on a sequence of examples. Our main motivation is the acceleration of the optimization of the regularized Cox partiallikelihood (the core model in survival analysis), but other settings can be considered as well.

The proposed algorithm is doubly stochastic in the sense that gradient steps are done using stochastic gradient descent (SGD) with variance reduction, and the inner expectations are approximated by a Monte-Carlo Markov-Chain (MCMC) algorithm. We derive conditions on the MCMC number of iterations guaranteeing convergence, and obtain a linear rate of convergence under strong convexity and a sublinear rate without this assumption.

We illustrate the fact that our algorithm improves the state-of-the-art solver for regularized Cox partial-likelihood on several datasets from survival analysis.

Introduction

During the past decade, advances in biomedical technology have brought high dimensional data to biostatistics and survival analysis in particular. Today's challenge for survival analysis lays in the analysis of massively high dimensional (numerous covariates) and large-scale (large number of observations) data, see in particular [START_REF] Murdoch | The inevitable application of big data to health care[END_REF]. Areas of application outside of biostatistics, such as economics (see [START_REF] Einav | Economics in the age of big data[END_REF]), or actuarial sciences (see [START_REF] Richards | A handbook of parametric survival models for actuarial use[END_REF]) are also concerned.

One of the core models of survival analysis is the Cox model (see [START_REF] David | Regression models and life tables (with discussion)[END_REF]) for which we propose, in the present paper, a novel scalable optimization algorithm tuned to handle II. Large-scale Cox model massively high dimensional and large-scale data. Survival data (y

i , x i , ± i ) n pat i =1 contains, for each individual i = 1, . . . , n pat , a features vector x i 2 R d , an observed time y i 2 R + , which is a failure time if ± i = 1 or a right-censoring time if ± i = 0. If D = {i : ± i = 1}
is the set of patients for which a failure time is observed, if n = |D| is the total number of failure times, and if

R i = { j : y j ∏ y i }
is the index of individuals still at risk at time y i , the negative Cox partial log-likelihood writes

°`(µ) = 1 n X i 2D h °x> i µ + log ≥ X j 2R i exp(x > j µ) ¥i (1)
for parameters µ 2 R d . This model can be regarded as a regression of the n failure times, using information from the n pat patients that took part to the study. With high-dimensional data, a regularization term is added to the partial likelihood to automatically favor sparsity in the estimates, see [T + 97] and [START_REF] Simon | Regularization paths for cox's proportional hazards model via coordinate descent[END_REF] for a presentation of Lasso and elastic-net penalizations, see also the review paper by [START_REF] Witten | Survival analysis with high-dimensional covariates[END_REF] for an exhaustive presentation. Several algorithms for the Cox model have been proposed to solve the regularized optimization problem at hand, see [PH07, SKJP09, Goe10] among others. These implementations use Newton-Raphson iterations, i.e. large matrices inversions, and can therefore not handle large-scale data. Cyclical coordinate descent algorithms have since been proposed and successfully implemented in R packages coxnet and fastcox, see [START_REF] Simon | Regularization paths for cox's proportional hazards model via coordinate descent[END_REF][START_REF] Yang | A cocktail algorithm for solving the elastic net penalized cox's regression in high dimensions[END_REF]. More recently [START_REF] Mittal | Large-scale parametric survival analysis[END_REF] adapted the column relaxation with logistic loss algorithm of [START_REF] Zhang | The value of unlabeled data for classiAEcation problems[END_REF] to the Cox model. The fact that all these algorithms are of cyclic coordinate descent type solve the problem, supported by Newton-Raphson type algorithms, of large matrices inversions. Yet another computationnally costly problem, speciAEc to the Cox model, has not been fully addressed: the presence of cumulative sums (over indices j 2 R i ) in the Cox partial likelihood. This problem was noticed in [START_REF] Mittal | Large-scale parametric survival analysis[END_REF], where a numerical workaround exploiting sparsity is proposed to reduce the computational cost. The cumulative sum prevents from successfully applying stochastic gradient algorithms, which are however known for their e ciency to handle large scale generalized linear models: see for instance SAG by [START_REF] Schmidt | Minimizing AEnite sums with the stochastic average gradient[END_REF], SAGA by [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], Prox-SVRG by [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] and SDCA by [START_REF] Shalev-Shwartz | Proximal stochastic dual coordinate ascent[END_REF] that propose very e cient stochastic gradient algorithms with constant step-size (hence achieving linear rates), see also Catalyst by [START_REF] Lin | A universal catalyst for AErst-order optimization[END_REF] that introduces a generic scheme to accelerate and analyze the convergence of those algorithms.

Such recent stochastic gradient algorithms have shown that it is possible to improve upon proximal full gradient algorithms for the minimization of convex problems of the form

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = 1 n n X i =1 f i (µ), (2)
where the functions f i are gradient-Lipschitz and h is prox-capable. These algorithms take advantage of the AEnite sum structure of f , by using some form of variance-reduced stochastic gradient descent. It leads to algorithms with a much smaller iteration complexity, as compared to proximal full gradient approach (FG), while preserving (or even improving) the linear convergence rate of FG in the strongly convex case. However, such algorithms are relevant 2. Comparison with previous work when gradients r f i have a numerical complexity much smaller than r f , such as for linear classiAEcation or regression problems, where r f i depends on a single inner product x > i µ between features x i and parameters µ. In this paper, motivated by the important example of the Cox partial likelihood (1), we consider the case where gradients r f i can have a complexity comparable to the one of r f . More precisely, we assume that they can be expressed as expectations, under a probability measure º i µ , of random variables G i (µ), i.e.,

r f i (µ) = E G i (µ)ªº i µ [G i (µ)]. (3)
This paper proposes a new doubly stochastic proximal gradient descent algorithm (2SVRG), that leads to a low iteration complexity, while preserving linear convergence under suitable conditions for problems of the form (2) + (3).

Our main motivation for considering this problem is to accelerate the training-time of the the penalized Cox partial-likelihood. The function °`(µ) is convex (as a sum of linear and log-sum-exp functions, see Chapter 3 of [START_REF] Boyd | Convex optimization[END_REF], and AEts in the setting (2) + (3). Indeed, AEx i 2 D and introduce

f i (µ) = °x> i µ + log ≥ X j 2R i exp(x > j µ) ¥ , so that r f i (µ) = °xi + X j 2R i x j º i µ ( j )
where

º i µ ( j ) = exp(x > j µ) P j 0 2R i exp(x > j 0 µ) , 8 j 2 R i .
This entails that r f i (µ) satisAEes (3) with G i (µ) a random variable valued in

{°x i + x j : j 2 R i }
and such that

P(G i (µ) = °xi + x j ) = º i µ ( j )
for j 2 R i . Note that the numerical complexity of r f i can be comparable to the one of r f , when y i is close to min

i y i (recalling that R i = { j : y j ∏ y i }).
Note also that a computational trick allows to compute r f (µ) with a complexity O(nd). Indeed, once all data points are sorted, the sum can be computed recursively. This makes this setting quite di erent from the usual case of empirical risk minimization (linear regression, logistic regression, etc.), where all the gradients r f i share the same low numerical cost.

Comparison with previous work SGD techniques. Recent proximal stochastic gradient descent algorithms by [DBLJ14], [XZ14], [SSZ12]

and [START_REF] Schmidt | Minimizing AEnite sums with the stochastic average gradient[END_REF] build on the idea of [START_REF] Robbins | A stochastic approximation method[END_REF] and [KW + 52]. Such algorithms are designed to tackle large-scale optimization problems (n is large), where it is assumed implicitly that the r f i (smooth gradients) have a low computational cost compared to r f , and where h is eventually non-di erentiable and is dealt with using a backward or projection step using its proximal operator.

The principle of SGD is, at each iteration t , to sample uniformly at random an index

i ª U [n]
, and to apply an update step of the form

µ t +1 √ µ t °∞t r f i (µ t ).
This step is based on an unbiased but very noisy estimate of the full gradient r f , so the choice of the step size ∞ t is crucial since it has to be decaying to curb the variance introduced by random sampling (excepted for averaged SGD in some particular cases, see [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n)[END_REF]). This tends to slow down convergence to a minimum µ ? 2 argmin µ2R d f (µ). Gradually reducing the variance of r f i for i ª U [n] as an approximation of r f allows to use larger -even constant -step sizes and to obtain faster convergence rates. This is the underlying idea of two recent methods -SAGA and SVRG respectively introduced in [DBLJ14], [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] -that use updates of the form

w t +1 √ µ t °∞≥ r f i (µ t ) °rf i ( μ) + 1 n n X j =1 r f j ( μ) ¥ , and
µ t +1 √ prox ∞h (w t +1
). In [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], μ is fully updated after a certain number of iterations, called phases, whereas in [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], μ is partially updated after each iteration. Both methods use stochastic gradient descent steps, with variance reduction obtained via the centered control variable °r

f i ( μ) + 1 n P n j =1 r f j ( μ), and achieve linear convergence when F is strongly- convex, namely EF (µ k ) °min 2R d F (x) = O(Ω k
) with Ω < 1, which make these algorithms state-of-the-art for many convex optimization problems. Some variants of SVRG [XZ14] also approximate the full gradient 1 n P n j =1 r f j ( μ) using mini-batchs to decrease the computing time of each phase, see [LJ17, HAV + 15].

Numerically hard gradients.

A very di erent, nevertheless classical, "trick" to reduce the complexity of the gradient computation, is to express it, whenever the statistical problem allows it, as the expectation, with respect to a non-uniform distribution º µ , of a random variable G(µ), i.e., r f

(µ) = E G(µ)ªº µ [G(µ)].
Optimization problems with such a gradient have generated an extensive literature from the AErst works by [START_REF] Robbins | A stochastic approximation method[END_REF], and [KW + 52]. Some algorithms are designed to construct stochastic approximations of the sub-gradient of f + h, see [NJLS09, JN + 11, Lan12, DHS11]. Others are based on proximal operators to better exploit the smoothness of f and the properties of h, see [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF][START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF][START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF]. In this paper, we shall focus on the second kind of algorithms. Indeed, our approach is closer to the one developed in [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF], though, as opposed to ours, the algorithm developed in this latter work is based on proximal full gradient algorithms (not doubly stochastic as ours) and does not guarantee a linear convergence.

Contrastive divergence. The idea to approximate the gradient using MCMC already appeared in the litterature of Undirected Graphical Models under the name of Contrastive Divergence, see [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF][START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF][START_REF] Carreira-Perpinan | On contrastive divergence learning[END_REF]. Indeeed, for this class of model, the gradient of the log-likelihood r f (µ) can be written as the di erence of two expectations: one -tractable 3. A doubly stochastic proximal gradient descent algorithm -with respect to the data discrete distribution X, the other -intractable -with respect to the model-dependent distribution p(•, µ). The idea of Contrastive Divergence relies in the approximation of the intractable expectation using MCMC, with few iterations of the chain. However, in the framework of Cox model, and also Conditional Random Fields (see Section 6 below), this is the gradient r f i (µ) that writes as an time-consuming expectation, see Equation 3.

Our setting. The setting of our paper is original in the sense that it combines both previous settings, namely stochastic gradient descent and MCMC. As in the stochastic gradient setting, the gradient can be expressed as the sum of n components, where n can be very large. However, since these components are time-consuming to compute directly, following the expectation based gradient computation setting, they are expressed as averaged values of some random variables. More precisely, the gradient r f i (µ) is replaced by an approximation b r f i (µ) obtained by an MCMC algorithm. Our algorithm is, to the best of our knowledge, the AErst one to propose a combination of two stochastic approximations in this way, hence the name doubly stochastic, which allow to deal with both, eventual large values for n and the inner complexity of each gradient r f i computation. The idea to mix SGD and MCMC has also been raised recently in the very di erent setting of implicit stochastic gradient descent, see [START_REF] Toulis | Implicit stochastic gradient descent for principled estimation with large datasets[END_REF]. Note also that in our approach we make two stochastic approximations to the gradient using random training points, while the doubly stochastic approach from [DXH + 14] performs two stochastic approximations to the gradient using random training points and random features for kernel methods.

A doubly stochastic proximal gradient descent algorithm

Our algorithm 2SVRG is built upon the algorithm SVRG via an approximation function ApproxMCMC. We AErst present the meta-algorithm without specifying the approximation function, and then provide two examples for ApproxMCMC.

2SVRG: a meta-algorithm

Following the ideas presented in the previous section, we design a doubly stochastic proximal gradient descent algorithm (2SVRG), by combining a variance reduction technique for SGD given by Prox-SVRG [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], and a Monte-Carlo Markov-Chain algorithm to obtain an approximation of the gradient r f j (µ) at each step. Thus, in the considered setting the full gradient writes

r f (µ) = E i ªU [r f i (µ)] = E i ªU E G i (µ)ªº i µ [G i (µ)],
where U is the uniform distribution on {1, . . . , n}, so our algorithm contains two levels of stochastic approximation: uniform sampling of i (the variance-reduced SGD part) for the AErst expectation, and an approximation of the second expectation w.r.t º i µ by means of Monte-Carlo simulation. The 2SVRG algorithm is described in Algorithm 6.

Following Prox-SVRG by [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], this algorithm decomposes in phases: iterations within a phase apply variance reduced stochastic gradient steps (with a backward proximal step, see Algorithm 6 Doubly stochastic proximal gradient descent (2SVRG)

1: Require: Number of phases K ∏ 1, phase-length m ∏ 1, step-size ∞ > 0, MCMC number of itera- tions per phase (N k ) K k=1 , starting point µ 0 2 R d 2: Initialize: μ √ µ 0 and compute r f i ( μ) for i = 1, . . . , n 3: for k = 1 to K do 4: for t = 0 to m °1 do 5: Pick i ª U [n] 6: b r f i (µ t ) √ ApproxMCMC(i , µ t , N k ) 7: d t = b r f i (µ t ) °rf i ( μ) + 1 n P n j =1 r f j ( μ) 8: ! t +1 √ µ t °∞d t 9: µ t +1 √ prox ∞h (! t +1 ) 10: end for 11: Update μ √ 1 m P m t =1 µ t , µ 0 √ μ, μk √ μ 12:
Compute r f i ( μ) for i = 1, . . . , n 13: end for 14: Return: μK lines 7 and 8 in Algorithm 6). At the end of a phase, a full-gradient is computed (lines 10, 11) and used in the next phase for variance reduction. Within a phase, each inner iteration samples uniformly at random an index i (line 4) and obtains an approximation of the gradient r f i at the previous iterate µ t by applying N k iterations of a Monte-Carlo Markov-Chain (MCMC) algorithm.

Intuitively, the sequence N k should be increasing with the phase number k, as we need more and more precision as the iterations goes on (this is conAErmed in Section 4). The important point of our algorithm resides precisely in this aspect: very noisy estimates can be used in the early phases of the algorithm, hence allowing for an overall low complexity as compared to a full gradient approach.

Choice of ApproxMCMC

We focus now on two implementations of the function ApproxMCMC based on two famous MCMC algorithms: Metropolis-Hastings and Importance Sampling.

Independent Metropolis-Hastings

When the º i µ are Gibbs probability measures, as for the previously described Cox partial loglikelihood (but for other models as well, such as Conditional Random Fields, see [LMP + 01]), one can apply Independent Metropolis-Hastings (IMH), see Algorithm 7 below, to obtain approximations b r f i of the gradients. In this case the produced chain is geometrically uniformly ergodic, see [START_REF] Robert | Monte carlo methods[END_REF], and therefore meets the general assumptions required in our results (see Proposition 1 below). The IMH algorithm uses a proposal distribution Q which is independent of the current state j l of the Markov chain. In the case of the Cox partial log-likelihood, at iteration t of phase k of Algorithm 6, we set º = º i µ t , and Q to be the uniform distribution over the set R i . We implemented two versions Algorithm 7 Independent Metropolis-Hastings (IMH) estimator (for the Cox model)

Require: Proposal distribution Q = U {R i }, starting point j 0 2 R i , stationary distribution º = º i µ t for l = 0, . . . , N k °1 do 1. Generate: j 0 ª Q. 2. Update: AE = min ≥ º( j 0 )Q( j l ) º( j l )Q( j 0
) ,1

¥ = min °exp((x j 0 °x j l ) > µ t
), 1 ¢ .

3. Take:

j l +1 = ( j 0 with probability AE j l otherwise. end for Return: °xi + 1 N k P N k l =1 x j l
of Algorithm 6 with IMH: one with a uniform proposal Q, the other one with an adaptative proposal e

Q. When we want to approximate r f i (µ), we can consider the adaptative proposal

e Q = º i μ,
where μ is the iterate we have computed at the end of the previous phase, see Line 10 of Algorithm 6. Since we compute the full gradient only once every phase, the probabilities º i μ( j ) are computed at the same time, which means that the use of an adaptative proposal adds no computational e ort. Morever, the theoretical guarantees given in Section 4 make no di erence between the two versions aformentionned, but a strong di erence is observed in practice.

Importance Sampling

To choice of the adaptative proposal above reduces the variance of the estimator given by Ap-proxMCMC. The idea of sampling with e Q = º i μ can also be used in an Importance Sampling estimator as well.

r f i (µ) = E G i (µ)ªº i µ [G i (µ)] = E G i (µ)ª e Q " G i (µ) º i µ (G i (µ)) e Q(G i (µ)) # Since the ratio º i µ (G i (µ))/ e Q(G i (µ))
still contains an expensive term to compute, we can divide the term above with

E e Q £ º i µ (G i (µ))/ e Q(G i (µ)) §
= 1 and approximate the resulting term. This trick provides an estimator called Normalized Importance Sampling estimator, which writes like this in the case of Cox partial likelihood:

b J N = N X k=1 (x j k °xi ) º i µ ( j k ) e Q( j k ) , N X k=1 º i µ ( j k ) e Q( j k ) , with j k ª e Q = °xi + N X k=1 exp((µ °μ) > x j k ) P N l =1 exp((µ °μ) > x j l ) x j k , with j k ª e Q
Section 4 below gives theoretical guarantees for Algorithm 6: linear convergence under strong-convexity of F is given in Theorem 1, and a convergence without strong convexity is given in Theorem 2. This improves the proximal stochastic gradient method of [AFM17],

II. Large-scale Cox model

Algorithm 8 Normalized Importance Sampling (NIS) estimator of r f i (µ) (for the Cox model)

Require: Proposal distribution e Q = º i μ, stationary distribution º i µ , V = 0 2 R d , S = 0 2 R for l = 1, . . . , N k do 1. Generate: j l ª e Q(•). 2. Update: V √ V + exp((µ °μ) > x j l )x j l . 3. Update: S √ S + exp((µ °μ) > x j l ).

end for Return

: °xi + V /S
where the best case rate is O(1/k 2 ) using Fista (see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]) acceleration scheme. Numerical illustrations are given in Section 5, where a fair comparison between several state-of-the-art algorithms is proposed.

Theoretical guarantees

DeAEnitions. All the functions f i and h are proper convex lower-semicontinuous on

R d . The norm k • k stands for the Euclidean norm on R d . A function f : R d ! R is L-smooth if it is di erentiable and if its gradient is L-Lipschitz, namely if kr f (x) °rf (y)k ∑ Lkx °yk for all x, y 2 R d . A function f : R d ! R is µ-strongly convex if f (x + y) ∏ f (x) + rf (x) > y + µ 2 kyk 2 for all x, y 2 R d i.e. if f °µ 2 k • k 2 is convex. The proximal operator of h : R d ! R is uniquely deAEned by prox h (x) = argmin y2R d {h(y) + 1 2 kx °yk 2 }.
Notations. We denote by i t the index randomly picked at the t th iteration, see line 4 in Algorithm 6. We introduce the error of the MCMC approximation

¥ t = b r f i t (µ t °1)°r f i t (µ t °1)
and the AEltration F

t = ae(µ 0 , i 1 , µ 1 ,...,i t , µ t
). In order to analyze the descent steps, we need di erent expectations: E t the expectation w.r.t the distribution of the pair (i

t , b r f i t (µ t °1)) con- ditioned on F t °1
, and E the expectation w.r.t all the random iterates (i

t , µ t
) of the algorithm. We also denote µ § = argmin µ2R d F (µ).

Assumptions. Assumption 1. We consider

F = f + h where f = 1 n P n i =1 f i , with each f i being convex and L i -smooth, L i > 0,
and h a lower semi-continuous and closed convex function. We denote L = max 1∑i ∑n L i . We assume that there exists B > 0 such that the iterates µ t satisfy sup

t ∏0 kµ t °µ § k ∑ B .
Assumption 2. We assume that the bias and the expected squared error of the Monte Carlo estimation can be bounded in the following way: Theorems. The theorems below provide upper bounds on the distance to the minimum in the strongly convex case, see Theorem 1 and in the convex case, see Theorem 2.

kE t ¥ t k ∑ C 1 N k and E t k¥ t k 2 ∑ C 2 N k ( 4
Theorem 1. Suppose that F = f + h is µ-strongly convex. Consider Algorithm 6, with a phase length m and a step-size ∞ 2 (0, 1 16L ) satisfying

Ω = 1 m∞µ(1 °8L∞) + 8L∞(1 + 1/m) 1 °8L∞ < 1.
(5)

Then, under Assumption 1 and Assumption 2, we have:

E[F ( μK )] °F (µ § ) ∑ Ω K ≥ F (µ 0 ) °F (µ § ) + K X l =1 D Ω l N l ¥ , (6)
where

D = 3∞C 2 +BC 1 1°8L∞ . In Theorem 1, the choice N k = k AE Ω °k with AE > 1 gives E[F ( μK )] °F (µ § ) ∑ D 0 Ω K where D 0 = F (µ 0
) °F (µ § ) + D P k∏1 k °AE and D > 0 is a numerical constant. This entails that 2SVRG achieves a linear rate under strong convexity.

Remark 1 (An important remark). The number N k of MCMC iterations is growing quickly with the phase number k. So, we use in practice an hybrid version of 2SVRG called HSVRG: 2SVRG is used for the AErst phases (usally 4 or 5 phases in our experiments), and as soon as N

k exceeds n, we switch to a mini-batch version of Prox-SVRG (SVRG-MB), see [START_REF] Nitanda | Stochastic proximal gradient descent with acceleration techniques[END_REF]. A precise description of HSVRG is given in Algorithm 9 from Section 5 below. Note that overall linear convergence of HSVRG is still guaranteed, since both 2SVRG and SVRG-MB decrease linearly the objective from one phase to the other. Theorem 2. Consider Algorithm 6, with a phase length m and a step-size ∞ 2 (0, 1 8L(2m+1) ). Then, under Assumption 1 and Assumption 2, we have:

E[F ( μK )] °F (µ § ) ∑ D 1 K + D 2 K K +1 X k=1 1 N k , (7)
where D 1 and D 2 depend on the constants of the problem, and where μK is the average of iterates μk until phase K .

In Theorem 2, the choice

N k = k AE with AE > 1 gives E[F ( μK )] °F (µ § ) ∑ D 3 K for a constant D 3 > 0.
This result is an improvement of the Stochastic Proximal Gradient algorithm from [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF] since it is not necessary to design a weighted averaged but just a simple average to reach the same convergence rate. Also, it provides a convergence guarantee for the non-strongly convex case, which is not proposed in [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF].

Theorems 1 and 2 show a trade-o between the linear convergence of the variancereduced stochastic gradient algorithm and the MCMC approximation error. The next proposition proves that Algorithm 7 satisAEes Assumption 2 under a general assumption on the proposal and the stationary distribution.

Proposition 1. Suppose that there exists M > 0 such that the proposal Q and the stationary distribution º satisfy º(x) ∑ MQ(x), for all x in the support of º. Then, the error ¥ t obtained by Algorithm 7 satisAEes Assumption 2.

Remark 2 (SpeciAEcs for the Cox partial likelihood). Note that the assumptions required in Proposition 1 are met for the Cox partial likelihood: in this case, a simple choice is

M = n max x2supp(º) º(x),
and the Monte Carlo error ¥ t induced by computing the gradient of f i at phase k using Algorithm 7 satisAEes (4) with

C 1 = 2 |R i | max j 2R i º i µ t °1 ( j ) C 2 = 36C 2 C 2 1 (1 +C 1 ) max j 2R i kx j k 2 2 ,
where C 2 is the Rosenthal constant of order 2, see Proposition 12 in [FM + 03].

Numerical experiments

We compare several solvers for the minimization of the objective given by an elastic-net penalization of the Cox partial likelihood

F (µ) = °`(µ) + ∏ ≥ AEkµk 1 + 1 °AE 2 kµk 2 2 ¥ ,
where we recall that the partial likelihood `is deAEned in Equation (1) and where ∏ > 0 and AE 2 [0, 1] are tuning parameters.

A fair comparison of algorithms. The doubly stochastic nature of the considered algorithms makes it hard to compare them to batch algorithms in terms of iteration number or epoch number (number of full passes over the data), as this is usually done for SGD-based algorithm. Hence, we proceed by plotting the evolution of F ( μ) °F (µ § ) (where

µ § 2 argmin u2R d F (u)
and μ is the current iterate of a solver) as a function of the number of inner products between a feature vector x i and µ, e ectively computed by each algorithm, to obtain the current iterate μ. This gives a fair way of comparing the e ective complexity of all algorithms. These algorithms however need a good starting point (near the actual minimizer) to achieve convergence (this fact is due to a diagonal approximation of the Hessian matrix, see [START_REF] Hastie | Generalized additive models[END_REF], Chapter 8.). They are therefore tuned to provide good path of solutions while varying by small steps the penalization parameter ∏. Indeed in this case, this starting point is naturally set at the minimizer at the previous value of ∏, when minimizing along a path but cannot be guessed outside of a path. We illustrate this fact on Hybrid SVRG algorithm Since N k exponentially increases, the 2SVRG's complexity is higher than SVRG's original complexity. However, the algorithm 2SVRG is very e cient during the AErst phases: we introduce an hybrid solver that begins with 2SVRG and switches to SVRG with mini-batchs (denoted SVRG-MB). Mini-batching simply consists in replacing single stochastic gradients r f i by an average over a subset B of size n mb uniformly selected at random. This is useful in our case, since we can use a computational trick (recurrence formula) to compute mini-batched gradients. In our experiments, we used n mb = 0.1n or n mb = 0.01n, a constant step-size ∞ designed for each dataset, and switched from 2SVRG to SVRG-MB after K S = 5 phases. We set

N k = n k/(K S +2) so that N k never exceeds n.
Baselines. We describe in this paragraph the algorithm that we put in competition in our experiments.

Algorithm 9 Hybrid SVRG (HSVRG)

1: Require: Number of phases before switching K S ∏ 1, total number of phases

K ∏ K S , phase-length m ∏ 1, step-size ∞ > 0, MCMC number of iterations per phase (N k ) K k=1 , starting point µ 0 2 R d 2: Initialize: μ √ µ 0 and compute r f i ( μ) for i = 1, . . . , n 3: for k = 1 to K S do 4:
for t = 0 to m °1 do 5:

Pick i ª U [n] 6: b r f i (µ t ) √ ApproxMCMC(i , µ t , N k ) 7: d t = b r f i (µ t ) °rf i ( μ) + 1 n P n j =1 r f j ( μ) 8: ! t +1 √ µ t °∞d t 9: µ t +1 √ prox ∞h (! t +1 ) 10: end for 11: Update μ √ 1 m P m t =1 µ t , µ 0 √ μ, μk √ μ 12: Compute r f i ( μ) for i = 1, . . . ,
for t = 0 to m mb °1 = b(m °1)/n mb c do 16: Pick a set of random indices B ª (U [n]) n mb 17: d t = rf B (µ t ) °rf B ( μ) + 1 n P n j =1 r f j ( μ) 18: ! t +1 √ µ t °∞d t 19: µ t +1 √ prox ∞h (! t +1 ) 20: end for 21: Update μ √ 1 m mb P m mb t =1 µ t , µ 0 √ μ, μk √ μ 22: end for 23: Return: μK
FISTA This is accelerated proximal gradient from [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] with backtracking linesearch. Inner products necessary inside the backtracking are counted as well.

L-BFGS-B

A state-of-the-art quasi-Newton solver which provides a usually strong baseline for many batch optimization algorithms, see [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF]. We use the original implementation of the algorithm proposed in python's scipy.optimize module. Nondi erentiability of the `1-norm in the elastic-net penalization is dealt with the standard trick of reformulating the problem, using the fact that |a| = a + + a °for a 2 R.

HSVRG-UNIF-IMH

Q = º • μ.
HSVRG-AIS This is Algorithm 9 where ApproxMCMC is done via Algorithm 8, that is Adaptative Importance Sampling.

SVRG-MB

Mini-Batch Prox-SVRG described in [START_REF] Nitanda | Stochastic proximal gradient descent with acceleration techniques[END_REF], which can be seen as Algorithm 9 (see below) with K S = 0. This is a simply stochastic algorithm, since there is no MCMC 6. Conclusion approximation of the gradients r f i . The question of mini-batch sizing is critical and is adressed in Section 10. We used n mb = 0.1n or n mb = 0.01n in our experiments.

The "simply stochastic" counterpart SVRG-MB is way slower than the corresponding doubly stochatic versions, since they rely on many computations of stochastic gradients r f i , which are numerically costly, as explained above. The same settings are used throughout all experiments, some of them being tuned by hand: steps size for the variants of HSVRG are taken as ∞ t = ∞ 0 2 {10 °2, 10 °3, 10 °4} where ∞ 0 depends on the dataset, the phase length m is equal to the number n of failures of each datasets as suggested in [START_REF] Konečn | Mini-batch semi-stochastic gradient descent in the proximal setting[END_REF]. As mentionned above, the doubly stochastic algorithms use di erent verions of ApproxMCMC.

Datasets

We compare algorithms on the following datasets. The AErst three are standard benchmarks in survival analysis, the fourth one is a large simulated dataset where the number of observations n exceeds the number of features d . This di ers from supervised gene expression data: such a large-scale setting happens for longitudinal clinical trials, medical adverse event monitoring and business data minings tasks. • We generated a Gaussian features matrix X with n = 10, 000 observations and d = 500 predictors, with a Toeplitz covariance and correlation equal to 0.5. The failure times follow a Weibull distribution. See Section 9 for details on simulation in this model.

We compare in Figures II.2 and II.3 all algorithms for ridge penalization, namely AE = 0 and ∏ = 1/ p n. Experiences with other values of AE and ∏ are given in Section (including the Lasso penalization for instance).

Conclusions.

The experiments AErst show that the solvers HSVRG-ADAP-IMH and HSVRG-AIS give better results than HSVRG-UNIF-IMF. However, the HSVRG solvers behave particularly well during the AErst phases where the gradients can be noisy -due to a small number of iterations of the MCMC -and still point a decent descent direction.

Conclusion

We have proposed a doubly stochastic gradient algorithm to extend SGD-like algorithms beyond the empirical risk minimization setting. The algorithm we proposed is the result of two di erent ideas: sampling from uniform distribution to avoid the computation of a large sum, and sampling using MCMC methods to avoid the computation of a more complicated expectation. We have also provided theoretical guarantees of convergence for both the convex and the strongly-convex setting.

This doubly stochastic gradient algorithm is very e cient during the early phases. The hybrid version of our algorithm, at the crossing of simply and doubly stochastic gradient algorithms, signiAEcantly outperforms state-of-the-art methods.

In a future work, we intend to extend our algorithm to Conditional Random Fields (CRF), where each subfunction's gradient takes the form

r f i (µ) = r(°log(p(y i |x i , µ)) = X Y 2Y i e H (X i ,Y ) > µ P Y 0 2Y i e H (X i ,Y 0 ) > µ (H (X i , Y ) °H (X i , Y i )),
for a certain function H (see Page 2 in [SBA + 15]). Notice that the Cox negative partial likelihood can be seen as a particular case of CRF by setting We AErst prove Proposition 1 that ensures that Algorithm 7 provides the bounds of Assumption 2.

X i = [x j ] j 2R i 2 R d £|R i | , Y i = [ j 2R i ] j 2R i 2 {0, 1} |R i | , H (X , Y ) = X Y and Y i = {[ j =k ] j 2R i : k 2 R i }.
Proof. Since there exists M > 0 such that the proposal Q and the stationary distribution º satisfy º(x) ∑ MQ(x), for all x in the support of º, the Theorem 7.8 in [START_REF] Robert | Monte carlo methods[END_REF] states that the Algorithm 7 produces a geometrically ergodic Markov kernel P with ergodicity constants uniformly controlled:

kP k (x, •) °ºk T V ∑ 2 µ 1 °1 M ∂ k , (8)
where P k is the kernel of the k th iteration of the algorithm and k • k T V is the total variation norm. Since b r f i t (µ t °1) is computed as the mean of the iterates of the Markov chain, a simple computation enables us to bound the bias of the error and Proposition 12 from [FM + 03] gives the upper bound for the expected squared error:

kE t ¥ t k ∑ C 1 N k and E t k¥ t k 2 ∑ C 2 N k (9)
where C 1 and C 2 are some AEnite constants, and N k the number of iterations of the Markov chain. It can be shown that C 1 = 2M and that C 2 is related to a constant from the Rosenthal's inequality.

Á

Preliminaries to the proofs of Theorems 1 and 2

In what follows, the key lemmas for the proofs of Theorems 1 and 2 are stated and proved when not directly borrowed from previous articles.

Lemma 1. For ¢ t := b r f i t (µ t °1) °rf i t ( μ) + rf ( μ) °rf (µ t °1
), we have:

E t k¢ t k 2 ∑ 8L[F (µ t °1) °F (µ § ) + F ( μ) °F (µ § )] + 3E t k¥ t k 2 .
The proof of Lemma 1 uses Lemma 1 in [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF].

Lemma 2. [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF][START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] Consider F satisfying Assumption 1. Then,

1 n n X i =1 kr f i (µ) °rf i (µ § )k 2 ∑ 2L[F (µ) °F (µ § )]
Proof of Lemma 1. For the sake of simplicity, we now denote

d t i = rf i (µ t °1) °rf i ( μ) and d t = r f (µ t °1) °rf ( μ), so that one gets ¢ t = d t i t °d t + ¥ t .
Then, using the expectation introduced in Section 4, we repeatedly use the identity

E t kªk 2 = E t kª °Et ªk 2 + kE t ªk 2 . First with ª = ¢ t (since E t d t i t = d t , one gets E t ª = E t ¥ t ) : E t k¢ t k 2 = E t kd t i t + ¥ t °°d t + E t ¥ t ¢ k 2 + kE t ¥ t k 2 then, successively with ª = d t i t + ¥ t , ª = d t + ¥ t and AEnally ª = ¥ t : E t k¢ t k 2 = E t kd t i t + ¥ t k 2 + kE t ¥ t k 2 °kd t + E t ¥ t k 2 = E t kd t i t + ¥ t k 2 + kE t ¥ t k 2 °°E t kd t + ¥ t k 2 °Et k¥ t °Et ¥ t k 2 ¢ = E t kd t i t + ¥ t k 2 + E t k¥ t k 2 °Et kd t + ¥ t k 2 .
Now we remark that E t kd t + ¥ t k 2 ∏ 0, and the identity ka + bk 2 ∑ 2kak 2 + 2kbk 2 gives the majoration

E t k¢ t k 2 ∑ 2E t kd t i t k 2 + 3E t k¥ t k 2 . Now rewriting d t i t = rf i t (µ t °1) °rf i t (µ § ) + rf i t (µ § ) °rf i t ( μ)
, the same identity leads to

E t k¢ t k 2 ∑ 4E t kr f i t (µ t °1) °rf i t (µ § )k 2 + 4E t kr f i t ( μ) °rf i t (µ § )k 2 + 3E t k¥ t k 2 .
The desired result follows applying twice Lemma 2.

Á

When F is µ-strongly convex, the next Lemma (Lemma 3 in [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF]) provides a key lower bound.

Lemma 3. [XZ14] Consider

F = f + h satifying Assumption 1, where f is L f -smooth, L f > 0, f is µ f -strongly convex, µ f ∏ 0, h is µ h -strongly convex, µ h ∏ 0. For any x, v 2 R d , we deAEne x + = prox ∞h (x °∞v), g = 1 ∞ (x °x+ ), where ∞ 2 (0, 1 L f ]. Then, for any y 2 R d : F (y) ∏ F (x + ) + g > (y °x) + ∞ 2 kg k 2 + µ f 2 ky °xk 2 + µ h 2 ky °x+ k 2 + (v °rf (x)) > (x + °y). (10)
Remark 3. Note that in Lemma 3, one can freely choose µ f and µ h (in particular one can take

µ f = 0 or µ h = 0), as long as µ f + µ h = µ.
The following Lemma comes from [AFM17] (Lemma 14):

Lemma 4. [AFM17] Consider F = f + h satifying Assumption 1, where f is L f -smooth, and T ∞ : x 7 ! prox ∞h [x °∞r f (x)] with ∞ 2 (0, 2/L f ]. Let x, y 2 R d , we have: kT ∞ (x) °T∞ (y)k ∑ kx °yk 7. Proofs

Proof of Theorem 1

Proof. The proof begins with the study of the distance kµ t °µ § k 2 between the phases k °1 and k. To ease the reading, when staying between these two phases, we write μ instead of μk°1 .

Introducing g t = 1 ∞ (µ t °1 °µt
), we may write:

kµ t °µ § k 2 = kµ t °1 °∞g t °µ § k 2 = kµ t °1 °µ § k 2 °2∞(g t ) > (µ t °1 °µ § ) + ∞ 2 kg t k 2 .
To upper bound the term °2∞(g t

) > (µ t °1 °µ §
)+∞ 2 kg t k 2 , we apply the Lemma 3 with x = µ t °1, x + = µ t and y = µ § . With again

¢ t = b r f i t (µ t °1) °rf i t ( μ) + rf ( μ) °rf (µ t °1), we obtain °(g t ) > (µ t °1 °µ § ) + ∞ 2 kg t k 2 ∑ F (µ § ) °F (µ t ) °µf 2 kµ t °1 °µ § k 2 °µh 2 kx t °µ § k 2 °(¢ t ) > (µ t °µ § ),
and

kµ t °µ § k 2 ∑ kµ t °1 °µ § k 2 + 2∞[F (µ § ) °F (µ t )] °2∞(¢ t ) > (µ t °µ § ). (11)
We now concentrate on the quantity °2∞(¢ t

) > (µ t °µ § ). Introducing ∫ t = prox ∞h [µ t °t 1 °∞r f (µ t °1)] 2 F t °1
i.e. the vector obtained from µ t °1 after an exact proximal gradient descent step, we get

°2∞(¢ t ) > (µ t °µ § ) = °2∞(¢ t ) > (µ t °∫t ) °2∞(¢ t ) > (∫ t °µ § ) ∑ 2∞k¢ t k • kµ t °∫t k °2∞(¢ t ) > (∫ t °µ § )
where the inequality follows from the Cauchy-Schwartz inequality. Now the non-expansiveness property of proximal operators kprox ∞h (x) °prox ∞h (y)k ∑ kx °yk leads to

°2∞(¢ t ) > (µ t °µ § ) ∑ 2∞k¢ t k • k{µ t °1 °∞(¢ t + rf (µ t °1))} °{µ t °1 °∞r f (µ t °1)}k °2∞(¢ t ) > (∫ t °µ § ) ∑ 2∞ 2 k¢ t k 2 °2∞(¢ t ) > (∫ t °µ § ).
Reminding that ∫ t 2 F t °1, we derive:

°2∞E t (¢ t ) > (µ t °µ § ) ∑ 2∞ 2 E t k¢ t k 2 °2∞(E t ¢ t ) > (∫ t °µ § ) ∑ 2∞ 2 E t k¢ t k 2 + 2∞kE t ¢ t k • k∫ t °µ § k,
the last inequality comes from the Cauchy-Schwartz inequality. Since µ § is the minimum of

F = f + h, it satisAEes µ § = prox ∞h [µ § °∞r f (µ § )]
. Thus, the Lemma 4 and the Assumption 1 on the sequence (µ t ) give us k∫ t °µ § k ∑ kµ t °1 °µ § k ∑ B . We also remark that E t ¢ t = E t ¥ t . For all t between phases k °1 and k, we AEnally apply Lemma 1 to obtain:

°2∞E t (¢ t ) > (µ t °µ § ) ∑ 16∞ 2 L[F (µ t °1) °F (µ § ) + F ( μ) °F (µ § )] + 6∞ 2 E t k¥ t k 2 + 2∞B kE t ¥ t k. ( 12
)
Taking the expectation E t on inequation (11) and combining with previous inequality leads to

E t kµ t °µ § k 2 ∑ kµ t °1 °µ § k 2 + 2∞[F (µ § ) °F (µ t )] + 16∞ 2 L[F (µ t °1) °F (µ § ) + F ( μ) °F (µ § )] + 6∞ 2 E t k¥ t k 2 + 2∞B kE t ¥ t k.
With the notation of Algorithm 6, μ = μk°1 = µ 0 . Now, applying iteratively the previous inequality over t = 1, 2, . . . , m and taking the expectation

E over i 1 , µ 1 , i 2 , µ 2 ,...,i
m , µ m , we obtain:

Ekµ m °µ § k 2 + 2∞[EF (µ m ) °F (µ § )] +2∞(1 °8L∞) m°1 X t =1 [EF (µ t ) °F (µ § )] ∑ kµ 0 °µ § k 2 + 16L∞ 2 [F (µ 0 ) °F (µ § ) + m(F ( μ) °F (µ § ))] + 6∞ 2 m X t =1 Ek¥ t k 2 + 2∞B m X t =1 EkE t ¥ t k.
Now, by convexity of F and the deAEnition μk = 1 m P m t =1 µ t , we may write F ( μk

) ∑ 1 m P m t =1 F (µ t
).

Noticing that 2∞(1 °8L∞) < 2∞ leads to

2∞(1 °8L∞)m[EF ( μk ) °F (µ § )] ∑ k μ °µ § k 2 + 16L∞ 2 (m + 1)[F ( μ) °F (µ § )] + 6∞ 2 m X t =1 Ek¥ t k 2 + 2∞B m X t =1 kE¥ t k.
Under the Assumption 2, we have

6∞ 2 m X t =1 Ek¥ t k 2 + 2∞B m X t =1 kE¥ t k ∑ (6∞ 2 C 2 + 2∞BC 1 ) m N k whereas the µ-strong convexity of F implies k μk°1 °µ § k 2 ∑ 2 µ [F ( μk°1 ) °F (µ §
)]. This leads to

EF ( μk ) °F (µ § ) ∑ Ω ≥ EF ( μk°1 ) °F (µ § ) ¥ + D N k
for D and Ω as deAEned in the theorem. Applying the last inequality recursively leads to the result.

Á 7. Proofs

Proof of Theorem 2

Proof. As at the begining of the proof of Theorem 1, we consider that we stand between phase k °1 and phase k of Algorithm 6 and consequently µ 0 = μk°1 . We use the same arguments until (11), with the di erence that, in this non-strongly convex case, we have µ

f = µ h = 0.
We obtain for all t between phases 1 and m

F (µ t ) °F (µ § ) ∑ 1 2∞ (kµ t °1 °µ § k 2 °kµ t °µ § k 2 ) °(µ t °µ § ) > ¢ t . Summing over t = 1, . . . , ø (for ø ∑ m) leads to ø X t =1 [F (µ t ) °F (µ § )] ∑ 1 2∞ ( ø°1 X t =0 kµ t °µ § k 2 °ø X t =1 kµ t °µ § k 2 ) °ø X t =1 (µ t °µ § ) > ¢ t . (13)
We now use Equation ( 13) (with ø = m) and the convexity of

k•k 2 with μk = 1 m P m t =1 µ t to write m X t =1 [F (µ t ) °F (µ § )] ∑ 1 2∞ √ m°1 X t =0 kµ t °µ § k 2 °mk μk °µ § k 2 ! °m X t =1 (µ t °µ § ) > ¢ t . (14)
Starting from Equation (13) again but now summing over l = 1, . . . , t , we get

1 2∞ (kµ 0 °µ § k 2 °kµ t °µ § k 2 ) °t X l =1 (µ l °µ § ) > ¢ l ∏ t X l =1 [F (µ l ) °F (µ § )] ∏ 0, (15)
where the last inequality follows from the deAEnition of µ § . In [START_REF]DXH +[END_REF], we now substitute kµ t °µ § k 2 by the upper bound derived from (15) to write (noticing that µ 0 = μk°1 ):

m X t =1 [F (µ t ) °F (µ § )] ∑ m 2∞ (k μk°1 °µ § k 2 °k μk °µ § k 2 ) °m°1 X t =1 t X l =1 (µ l °µ § ) > ¢ l °m X t =1 (µ t °µ § ) > ¢ t ∑ m 2∞ (k μk°1 °µ § k 2 °k μk °µ § k 2 ) °m X t =1 (m + 1 °t )(µ t °µ § ) > ¢ t .
As in the proof of Theorem 1 (see Equation ( 12)), each term °Et

(µ t °µ § ) > ¢ t is upper bounded by 8∞L[F (µ t °1) °F (µ § ) + F ( μk°1 ) °F (µ § )] + 3∞E t ||¥ t || 2 + B ||E t ¥ t ||. Now with m + 1 °t ∑ m and Assumption 2, we obtain: 1 m m X t =1 E[F (µ t ) °F (µ § )] ∑ 1 2∞ (k μk°1 °µ § k 2 °Ek μk °µ § k 2 ) + 8L∞ © m X t =1 [EF (µ t °1) °F (µ § )] + F (µ § ) °E[F (µ m )] + (m + 1)[E[F ( μk°1 )] °F (µ § )] ™ + m 3∞C 2 + BC 1 N k .
By deAEnition of ∞, we have 8Lm∞ < 1, and we can use the convexity of F to lower bound the left hand side. With the inequality

E[F (µ m )] °F (µ §
) ∏ 0, one has:

(1 °8L∞m) h E[F ( μk )] °F (µ § ) i ∑ 1 2∞ ≥ k μk°1 °µ § k 2 °Ek μk °µ § k 2 ¥ + 8L∞(m + 1) h E[F ( μk°1 )] °F (µ § ) i + m 3∞C 2 + BC 1 N k
We now take the expectation E on all iterates of the algorithm i.e. on the iterates i

1 , µ 1 , i 2 , µ 2 ,...,i m , µ m from the AErst phase. Introduce the notations A k = E[F ( μk )] °F (µ §
) and a = (8L∞(m + 1))/(1 °50 8. Supplementary experiments 8Lm∞) < 1 , last inequality leads to:

A k °a A k°1 ∑ 1 2∞(1 °8Lm∞) ≥ Ek μk°1 °µ § k 2 °Ek μk °µ § k 2 ¥ + D N k ,
where D is deAEned in the theorem. Summing over the phases k = 1, 2, . . . , K + 1 and lower bounding A K +1 with 0, we obtain:

(1 °a) K X k=1 A k ∑ a A 0 + 1 2∞(1 °8Lm∞) k μ0 °µ § k 2 + K +1 X k=1 D N k
The last argument is the use of the convexity of F . Remark the explicit forms of the constants in the theorem:

D 1 = a 1 °a A 0 + 1 1 °a k μ0 °µ § k 2 2∞(1 °8Lm∞) and D 2 = D 1°a .
Á

Supplementary experiments

We have tested all algorithms with other settings for the penalization. Namely, we considered:

High lasso. We take AE = 1 and ∏ = 1/

Simulation of data

With Cox model, the hazard ratio for the failure time T i of the i th patient takes the form:

∏ i (t ) = ∏ 0 (t ) exp(x > i µ),
where ∏ 0 (t ) is a baseline hazard ratio, and x i 2 R d the covariates of the i th patient. We AErst simulate the feature matrix X 2 R n£d as a Gaussian vector with a Toepliz covariance, where the correlation between features j and j 0 is equal to Ω | j °j 0 | , for some Ω 2 (0, 1). We want now to simulate the observed time y i that corresponds to x i . We denote the cumulative hazard function

§(t ) = R t 0 ∏(s)d s. Using the deAEnition ∏(t ) = f (t ) 1°F (t ) , we know that §(t ) = °log(1 °F (t ))
, where f is the p.d.f. and F is the c.d.f. of T . It is easily seen that §(T ) has distribution Exp(1) (Exponential with intensity equal to 1): since § is an increasing function, we have

P( §(T ) ∏ t ) = P(T ∏ § °1(t )) = Z 1 § °1(t ) f (s)d s = 1 °F ( § °1(t )) = exp(° §( § °1(t )))
= exp(°t ), so that simulating failure times is simply achieved by using

T i = § °1(E i ) where E i ª Exp(1).
To compute §, we should have a parametric form for ∏ 0 . We assume that T follows the Weibull distribution W (1, ∫) (when x i = 0). This choice is motivated by the following facts: • Its cumulative hazard function is easy to invert. Indeed the hazard ratio is given by

∏ 0 (t ) = ∫t ∫°1 e °t ∫ 1°(1°e °t ∫ ) = ∫t ∫°1 , so that § °1(y ) = ≥ y exp(x > i µ) ¥ 1/∫ .
• It enables two di erent trends -increasing or decreasing -for the baseline hazard ratio that correspond to two typical behaviours in the medical AEeld.

-decreasing: after taking a treatment, time before a side-e ect's appearence

u i ∑ k ∂ = n mb Y i =1 P(u i ∑ k) = µ bkc n ∂ n mb , P µ max i 2B |R i | ∑ ck ∂ = µ bkc n ∂ n mb , P µ max i 2B |R i | ∏ a ∂ = 1 °µ ba/cc n ∂ n mb , for a < n pat
The third equation leads us to consider 1 ø n mb ø n to prevent both max

Uncover Hawkes causality without parametrization 1 Introduction

In many applications, one needs to deal with data containing a very large number of irregular timestamped events that are recorded in continuous time. These events can reØect, for instance, the activity of users on a social network, see [SAD + 16], the high-frequency variations of signals in AEnance, see [START_REF] Bacry | Hawkes processes in AEnance[END_REF], the earthquakes and aftershocks in geophysics, see [START_REF] Ogata | Space-time point-process models for earthquake occurrences[END_REF], the crime activity, see [MSB + 11] or the position of genes in genomics, see [START_REF] Reynaud-Bouret | Adaptive estimation for hawkes processes; application to genome analysis[END_REF]. The succession of the precise timestamps carries a great deal of information about the dynamics of the underlying systems. In this context, multidimensional counting processes based models play a paramount role. Within this framework, an important task is to recover the mutual inØuence of the nodes (i.e., the di erent components of the counting process), by leveraging on their timestamp patterns, see, for instance, [BM16, LV14, LM11, ZZS13, GRLS13, FWR + 15, XFZ16].

Consider a set of nodes I = {1, . . . , d }. For each i 2 I , we observe a set Z i of events, where each ø 2 Z i labels the occurrence time of an event related to the activity of i . The events of all nodes can be represented as a vector of counting processes

N t = [N 1 t ••• N d t ] >
, where N i t counts the number of events of node i until time t 2 R + , namely

N i t = P ø2Z i {ø∑t } . The vector of stochastic intensities ∏ t = [∏ 1 t •••∏ d t ]
> associated with the multivariate counting process N t is deAEned as

∏ i t = lim d t!0 P(N i t +d t °N i t = 1|F t ) d t
for i 2 I , where the AEltration F t encodes the information available up to time t . The coordinate ∏ i t gives the expected instantaneous rate of event occurrence at time t for node i . The vector ∏ t characterizes the distribution of N t , see [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF], and patterns in the events time-series can be captured by structuring these intensities.

The Hawkes process introduced in [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF] corresponds to an autoregressive structure of the intensities in order to capture self-excitation and cross-excitation of nodes, which is a phenomenon typically observed, for instance, in social networks, see for instance [START_REF] Crane | Robust dynamic classes revealed by measuring the response function of a social system[END_REF]. Namely, N t is called a Hawkes point process if the stochastic intensities can be written as

∏ i t = µ i + d X j =1 Z t 0 ¡ i j (t °t 0 )d N j t 0 ,
where µ i 2 R + is an exogenous intensity and ¡ i j are positive, integrable and causal (with support in R + ) functions called kernels encoding the impact of an action by node j on the activity of node i . Note that when all kernels are zero, the process is a simple homogeneous multivariate Poisson process.

Most of the litterature uses a parametric approach for estimating the kernels. With no doubt, the most popular parametrization form is the exponential kernel ¡ i j

(t ) = AE i j Ø i j e °Øi j t
because it deAEnitely simpliAEes the inference algorithm (e.g., the complexity needed for computing the likelihood is much smaller). When d is large, in order to reduce the number of parameters, some authors choose to arbitrarily share the kernel shapes across the di erent nodes. Thus, for instance, in [YZ13, ZZS13, FWR + 15], they choose ¡ i j (t ) = AE i j h(t ) with AE i j 2 R + quantiAEes the intensity of the inØuence of j on i and h(t ) a (normalized) function that characterizes the time-proAEle of this inØuence and that is shared by all couples of nodes (i , j ) (most often, it is chosen to be either exponential h(t ) = Øe °Øt or power law h(t ) = Øt °(Ø+1) ). Both approaches are, most of the time, highly non-realistic. On the one hand there is a priori no reason for assuming that the time-proAEle of the inØuence of a node j on a node i does not depend on the pair (i , j ). On the other hand, assuming an exponential shape or a power law shape for a kernel arbitrarily imposes an event impact that is always instantly maximal and that can only decrease with time, while in practice, there may exist a latency between an event and its maximal impact.

In order to have more Øexibility on the shape of the kernels, nonparametric estimation can be considered. Expectation-Maximization algorithms can be found in [START_REF] Lewis | A nonparametric em algorithm for multiscale hawkes processes[END_REF] (for d = 1) or in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF] (d > 1). An alternative method is proposed in [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF] where the nonparametric 1. Introduction estimation is formulated as a numerical solving of a Wiener-Hopf equation. Another nonparametric strategy considers a decomposition of kernels on a dictionary of function h 1 ,...,h

K , namely ¡ i j (t ) = P K k=1 a i j k h k (t )
, where the coe cients a i j k are estimated, see [HRBR + 15, LV14] and [START_REF] Xu | Learning granger causality for hawkes processes[END_REF], where group-lasso is used to induce a sparsity pattern on the coe cients a i j k that is shared across k = 1, . . . , K .

Such methods are heavy when d is large, since they rely on likelihood maximization or least squares minimization within an over-parametrized space in order to gain Øexibility on the shape of the kernels. This is problematic, since the original motivation for the use of Hawkes processes is to estimate the inØuence and causality of nodes, the knowledge of the full parametrization of the model being of little interest for causality purpose.

Our paper solves this problem with a di erent and more direct approach. Instead of trying to estimate the kernels ¡ i j , we focus on the direct estimation of their integrals. Namely, we want to estimate the matrix G = [g i j ] where

g i j = Z +1 0 ¡ i j (u) du ∏ 0 for 1 ∑ i , j ∑ d . (1)
As it can be seen from the cluster representation of Hawkes processes ([HO74]), this integral represents the mean total number of events of type i directly triggered by an event of type j , and then encodes a notion of causality. Actually, as detailed below (see Section 2.1), such integral can be related to the Granger causality ( [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF]).

The main idea of the method we developed in this paper is to estimate the matrix G directly using a matching cumulants (or moments) method. Apart from the mean, we shall use second and third-order cumulants which correspond respectively to centered second and third-order moments. We AErst compute an estimation c M of these centered moments M (G) (they are uniquely deAEned by G). Then, we look for a matrix b G that minimizes the L 2 error kM ( b G) °c Mk 2 . Thus the integral matrix b G is directly estimated without making hardly any assumptions on the shape the involved kernels. As it will be shown, this approach turns out to be particularly robust to the kernel shapes, which is not the case of all previous Hawkes-based approaches that aim causality recovery. We call this method NPHC (Non Parametric Hawkes Cumulant), since our approach is of nonparametric nature. We provide a theoretical analysis that proves the consistency of the NPHC estimator. Our proof is based on ideas from the theory of Generalized Method of Moments (GMM) but requires an original technical trick since our setting strongly departs from the standard parametric statistics with i.i.d observations. Note that moment and cumulant matching techniques proved particularly powerful for latent topic models, in particular Latent Dirichlet Allocation, see [START_REF] Podosinnikova | Rethinking lda: moment matching for discrete ica[END_REF]. A small set of previous works, namely [START_REF] Fonseca | Hawkes process: Fast calibration, application to trade clustering, and di usive limit[END_REF][START_REF] Aït-Sahalia | Modeling AEnancial contagion using mutually exciting jump processes[END_REF], already used method of moments with Hawkes processes, but only in a parametric setting. Our work is the AErst to consider such an approach for a nonparametric counting processes framework.

The paper is organized as follows: in Section 2, we provide the background on the integrated kernels and the integrated cumulants of the Hawkes process. We then introduce the method, investigate its complexity and explain the consistency result we prove. In Section 3, we estimate the matrix of Hawkes kernels' integrals for various simulated datasets and for real datasets, namely the MemeTracker database and AEnancial order book data. We then provide in Section 4 the technical details skipped in the previous parts and the proof of our consistency result. Section 5 contains concluding remarks.

NPHC: The Non Parametric Hawkes Cumulant method

In this Section, we provide the background on integrals of Hawkes kernels and integrals of Hawkes cumulants. We then explain how the NPHC method enables estimating G.

Branching structure and Granger causality

From the deAEnition of Hawkes process as a Poisson cluster process, see [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] or [START_REF] Hawkes | A cluster process representation of a selfexciting process[END_REF], g i j can be simply interpreted as the average total number of events of node i whose direct ancestor is a given event of node j (by direct we mean that interactions mediated by any other intermediate event are not counted). In that respect, G not only describes the mutual inØuences between nodes, but it also quantiAEes their direct causal relationships. Namely, introducing the counting function N i √ j t that counts the number of events of i whose direct ancestor is an event of j , we know from [START_REF] Bacry | Hawkes processes in AEnance[END_REF] that

E[d N i √ j t ] = g i j E[d N j t ] = g i j § j d t, (2)
where we introduced § i as the intensity expectation, namely satisfying E[d N i t ] = § i d t. Note that § i does not depend on time by stationarity of N t , which is known to hold under the stability condition kGk < 1, where kGk stands for the spectral norm of G. In particular, this condition implies the non-singularity of I d °G.

Since the question of a real causality is too complex in general, most econometricians agreed on the simpler deAEnition of Granger causality [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF]. Its mathematical formulation is a statistical hypothesis test: X causes Y in the sense of Granger (u) = 0 for u 2 R + . Since the kernels take positive values, the latter condition is equivalent to

R 1 0 ¡ i j
(u)du = 0. In the following, we'll refer to learning the kernels' integrals as uncovering causality since each integral encodes the notion of Granger causality, and is also linked to the number of events directly caused from a node to another node, as described above at Eq. (2).

Integrated cumulants of the Hawkes process

A general formula for the integral of the cumulants of a multivariate Hawkes process is provided in [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF]. As explained below, for the purpose of our method, we only need to consider cumulants up to the third order. Given 1 ∑ i , j , k ∑ d , the AErst three integrated 2. NPHC: The Non Parametric Hawkes Cumulant method cumulants of the Hawkes process can be deAEned as follows thanks to stationarity:

§ i d t = E(d N i t )
(3)

C i j d t = Z ø2R ≥ E(d N i t d N j t +ø ) °E(d N i t )E(d N j t +ø ) ¥ (4) K i jk d t = ZZ ø,ø 0 2R 2 ≥ E(d N i t d N j t +ø d N k t +ø 0 ) + 2E(d N i t )E(d N j t +ø )E(d N k t +ø 0 ) °E(d N i t d N j t +ø )E(d N k t +ø 0 ) °E(d N i t d N k t +ø 0 )E(d N j t +ø ) °E(d N j t +ø d N k t +ø 0 )E(d N i t ) ¥ , (5)
where Eq. ( 3) is the mean intensity of the Hawkes process, the second-order cumulant (4) refers to the integrated covariance density matrix and the third-order cumulant (5) measures the skewness of N t . Using the martingale representation from [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF] or the Poisson cluster process representation from [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF], one can obtain an explicit relationship between these integrated cumulants and the matrix G. If one sets

R = (I d °G) °1, (6)
straightforward computations (see Section 4) lead to the following identities:

§ i = d X m=1 R i m µ m (7) C i j = d X m=1 § m R i m R j m (8) K i jk = d X m=1 (R i m R j m C km + R i m C j m R km +C i m R j m R km °2 § m R i m R j m R km ). (9)
Equations ( 8) and (9) are proved in Section 4. Our strategy is to use a convenient subset of Eqs.

(3), (4) and (5) to deAEne M, while we use Eqs. ( 7), (8) and (9) in order to construct the operator that maps a candidate matrix R to the corresponding cumulants M (R). By looking for b R that minimizes R 7 ! kM (R) °c Mk 2 , we obtain, as illustrated below, good recovery of the ground truth matrix G using Equation (6).

The simplest case d = 1 has been considered in [START_REF] Hardiman | Branching-ratio approximation for the self-exciting hawkes process[END_REF], where it is shown that one can choose M = {C 11 } in order to compute the kernel integral. Eq. ( 8) then reduces to a simple second-order equation that has a unique solution in R (and consequently a unique G) that accounts for the stability condition (kGk < 1).

Unfortunately, for d > 1, the choice M = {C i j } 1∑i ∑ j ∑d is not su cient to uniquely determine the kernels integrals. In fact, the integrated covariance matrix provides d (d + 1)/2 independent coe cients, while d 2 parameters are needed. It is straightforward to show that the remaining d (d °1)/2 conditions can be encoded in an orthogonal matrix O, reØecting the fact that Eq. ( 8) is invariant under the change R ! OR, so that the system is under-determined.

Our approach relies on using the third order cumulant tensor K = [K i jk ] which contains (d 3 + 3d 2 + 2d )/6 > d 2 independent coe cients that are su cient to uniquely AEx the matrix G. This can be justiAEed intuitively as follows: while the integrated covariance only contains symmetric information, and is thus unable to provide causal information, the skewness given by the third order cumulant in the estimation procedure can break the symmetry between past and future so as to uniquely AEx G. Thus, our algorithm consists of selecting d 2 third-order cumulant components, namely M = {K i i j } 1∑i , j ∑d . In particular, we deAEne the estimator of R as b R 2 argmin R L (R), where

L (R) = (1 °∑)kK c (R) °c K c k 2 2 + ∑kC (R) °b C k 2 2 , (10)
where k • k 2 stands for the Frobenius norm, K c = {K i i j } 1∑i , j ∑d is the matrix obtained by the contraction of the tensor K to d 2 indices, C is the covariance matrix, while c K c and b C are their respective estimators, see Equations ( 12), (13) below. It is noteworthy that the above mean square error approach can be seen as a particular Generalized Method of Moments (GMM), see [START_REF] Hall | Generalized method of moments[END_REF]. This framework allows us to determine the optimal weighting matrix involved in the loss function. However, this approach is unusable in practice, since the associated complexity is too high. Indeed, since we have d 2 parameters, this matrix has d 4 coe cients and GMM calls for computing its inverse leading to a O(d 6

) complexity. In this work, we use the coe cient ∑ to scale the two terms, as

∑ = k c K c k 2 2 k c K c k 2 2 + k b C k 2 2 ,
see Section 4.4 for an explanation about the link between ∑ and the weighting matrix. Finally, the estimator of G is straightforwardly obtained as

b G = I d °b R °1,
from the inversion of Eq. (6). Let us mention an important point: the matrix inversion in the previous formula is not the bottleneck of the algorithm. Indeed, its has a complexity O(d 3

)

that is cheap compared to the computation of the cumulants when n = max i |Z i | ¿ d , which is the typical scaling satisAEed in applications. Solving the considered problem on a larger scale, say d ¿ 10 3 , is an open question, even with state-of-the-art parametric and nonparametric approaches, see for instance [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF][START_REF] Xu | Learning granger causality for hawkes processes[END_REF][START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF][START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF], where the number of components d in experiments is always around 100 or smaller. Note that, actually, our approach leads to a much faster algorithm than the considered state-of-the-art baselines, see Tables 1-4 from Section 3 below.

Estimation of the integrated cumulants

In this section we present explicit formulas to estimate the three moment-based quantities listed in the previous section, namely, §, C and K . We AErst assume there exists H > 0 such that the truncation from (°1, +1) to [°H , H ] of the domain of integration of the quantities appearing in Eqs. (4) and (5), introduces only a small error. In practice, this amounts to neglecting border e ects in the covariance density and in the skewness density that is a good 2. NPHC: The Non Parametric Hawkes Cumulant method approximation if the support of the kernel ¡ i j (t ) is smaller than H and the spectral norm kGk satisAEes kGk < 1. In this case, given a realization of a stationary Hawkes process {N t : t 2 [0, T ]}, as shown in Section 4, we can write the estimators of the AErst three cumulants (3), (4) and (5) as

b § i = 1 T X ø2Z i 1 = N i T T (11) b C i j = 1 T X ø2Z i ≥ N j ø+H °N j ø°H °2H b § j ¥ (12) b K i jk = 1 T X ø2Z i ≥ N j ø+H °N j ø°H °2H b § j ¥ • ≥ N k ø+H °N k ø°H °2H b § k ¥ °b § i T X ø2Z j X ø 0 2Z k (2H °|ø 0 °ø|) + + 4H 2 b § i b § j b § k . ( 13
)
Let us mention the following facts.

Bias. While the AErst cumulant §i is an unbiased estimator of § i , the other estimators b C i j and b K i jk introduce a bias. However, as we will show, in practice this bias is small and hardly a ects numerical estimations (see Section 3). This is conAErmed by our theoretical analysis, which proves that if H does not grow too fast compared to T , then these estimated cumulants are consistent estimators of the theoretical cumulants (see Section 2.6).

Complexity.

The computations of all the estimators of the AErst, second and third-order cumulants have complexity respectively O(nd), O(nd 2 ) and O(nd 3 ), where n = max

i |Z i |.
However, our algorithm requires a lot less than that: it computes only d 2 third-order terms, of the form b K i i j , leaving us with only O(nd 2 ) operations to perform.

Symmetry. While the values of § i

,C i j and K i jk are symmetric under permutation of the indices, their estimators are generally not symmetric. We have thus chosen to symmetrize the estimators by averaging their values over permutations of the indices. Worst case is for the estimator of K c , which involves only an extra factor of 2 in the complexity.

The NPHC algorithm

The objective to minimize in Equation ( 10) is non-convex. More precisely, the loss function is a polynomial of R of degree 6. However, the expectations of cumulants § and C deAEned in Eq. ( 4) and (5) that appear in the deAEnition of L (R) are unknown and should be replaced with b § and b C . We denote f L (R) the objective function, where the expectations of cumulants § i and C i j have been replaced with their estimators in the right-hand side of Eqs. ( 8) and (9):

f L (R) = (1 °∑)kR Ø2 b C > + 2[R Ø ( b C °R b L)]R > °c K c k 2 2 + ∑kR b LR > °b C k 2 2 ( 14
)
As explained in [CHM + 15], the loss function of a typical multilayer neural network with simple nonlinearities can be expressed as a polynomial function of the weights in the network, whose degree is the number of layers. Since the loss function of NPHC writes as a polynomial of degree 6, we expect good results using optimization methods designed to train deep multilayer neural networks. We used the AdaGrad from [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF], a variant of the Stochastic Gradient Descent with adaptive learning rates. AdaGrad scales the learning rates coordinate-wise using the online variance of the previous gradients, in order to incorporate second-order information during training. The NPHC method is summarized schematically in Algorithm 10. Our problem being non-convex, the choice of the starting point has a major e ect on the convergence. Here, the key is to notice that the matrices R that match Equation ( 8

Complexity of the algorithm

Compared with existing state-of-the-art methods to estimate the kernel functions, e.g., the ordinary di erential equations-based (ODE) algorithm in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], the Granger Causalitybased algorithm in [START_REF] Xu | Learning granger causality for hawkes processes[END_REF], the ADM4 algorithm in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], and the Wiener-Hopf-based algorithm in [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF], our method has a very competitive complexity. This can be understood by the fact that those methods estimate the kernel functions, while in NPHC we only estimate their integrals. The ODE-based algorithm is an EM algorithm that parametrizes the kernel function with M basis functions, each being discretized to L points. The basis functions are updated after solving M Euler-Lagrange equations. If n denotes the maximum number of events per component (i.e. n = max 1∑i ∑d |Z i |) then the complexity of one iteration of the algorithm is O(Mn 3

d 2 + ML(nd + n 2
)). The Granger Causality-based algorithm is similar to the previous one, without the update of the basis functions, that are Gaussian kernels. The complexity per iteration is O(Mn 3 d 2

). The algorithm ADM4 is similar to the two algorithms above, as EM algorithm as well, with only one exponential kernel as basis function. The complexity per iteration is then O(n 3 d 2

). The Wiener-Hopf-based algorithm is not iterative, 2. NPHC: The Non Parametric Hawkes Cumulant method on the contrary to the previous ones. It AErst computes the empirical conditional laws on many points, and then invert the Wiener-Hopf system, leading to a O(nd 2

L + d 4 L 3
) computation. Similarly, our method AErst computes the integrated cumulants, then minimize the objective function with N iter iterations, and invert the resulting matrix b R to obtain b G. In the end, the complexity of the NPHC method is O(nd 2 + N iter d 3

). According to this analysis, summarized in Table III.1 below, one can see that in the regime n ¿ d , the NPHC method outperforms all the other ones.

O(N iter M (n 3 d 2 + L(nd + n 2 ))) GC [XFZ16] O(N iter Mn 3 d 2 ) ADM4 [ZZS13] O(N iter n 3 d 2 ) WH [BM16] O(nd 2 L + d 4 L 3 ) NPHC O(nd 2 + N iter d 3 )

Theoretical guarantee: consistency

The NPHC method can be phrased using the framework of the Generalized Method of Moments (GMM). GMM is a generic method for estimating parameters in statistical models. In order to apply GMM, we have to AEnd a vector-valued function g (X , µ) of the data, where X is distributed with respect to a distribution P µ 0 , which satisAEes the moment condition: E[g (X , µ)] = 0 if and only if µ = µ 0 , where µ 0 is the "ground truth" value of the parameter. Based on i.i.d. observed copies x 1 ,..., x n of X , the GMM method minimizes the norm of the empirical mean over n samples, k 1 n P n i =1 g (x i , µ)k, as a function of µ, to obtain an estimate of µ 0 .

In the theoretical analysis of NPHC, we use ideas from the consistency proof of the GMM, but the proof actually relies on very di erent arguments. Indeed, the integrated cumulants estimators used in NPHC are not unbiased, as the theory of GMM requires, but asymptotically unbiased. Moreover, the setting considered here, where data consists of a single realization {N t } of a Hawkes process strongly departs from the standard i.i.d setting. Our approach is therefore based on the GMM idea but the proof is actually not using the theory of GMM.

In the following, we use the subscript T to refer to quantities that only depend on the process (N t ) in the interval [0, T ] (e.g., the truncation term H T , the estimated integrated covariance b C T or the estimated kernel norm matrix b G T ). In the next equation, Ø stands for the Hadamard product and Ø2 stands for the entrywise square of a matrix. We denote

G 0 = I d °R°1
0 the true value of G, and the R 2d £d valued vector functions Theorem 1 (Consistency of NPHC). Suppose that (N t ) is observed on R + and assume that

g 0 (R) = ∑ C °RLR > K c °RØ2 C > °2[R Ø (C °RL)]R > ∏ b g T (R) = " b C T °R b L T R > c K c T °RØ2 b C > T °2[R Ø ( b C T °R b L T )]R > . # Using these notations, f L T (R)
1. g 0 (R) = 0 if and only if R = R 0 ; 2. R 2 £,
where £ is a compact set;

the spectral radius of the kernel norm matrix satisAEes kG

0 k < 1; 4. H T ! 1 and H 2 T /T ! 0. Then b G T = I d °µarg min R2£ f L T (R) ∂ °1 P ! G 0 .
The proof of the Theorem is given in Section 4.5 below. Assumption 3 is mandatory for stability of the Hawkes process, and Assumptions 3 and 4 are su cient to prove that the estimators of the integrated cumulants deAEned in Equations ( 11), ( 12) and (13) are asymptotically consistent. Assumption 2 is a very mild standard technical assumption allowing to prove consistency for estimators based on moments. Assumption 1 is a standard asymptotic moment condition, that allows to identify parameters from the integrated cumulants.

Numerical Experiments

In this Section, we provide a comparison of NPHC with the state-of-the art, on simulated datasets with di erent kernel shapes, the MemeTracker dataset (social networks) and the order book dynamics dataset (AEnance).

Simulated datasets.

We simulated several datasets with Ogata's Thinning algorithm [START_REF] Ogata | On lewis' simulation method for point processes[END_REF] using the open-source library tick1 , each corresponding to a shape of kernel: rectangular, exponential or power law kernel, see Figure III.1 below.

The integral of each kernel on its support equals AE, 1/Ø can be regarded as a characteristic time-scale and ∞ is the scaling exponent for the power law distribution and a delay parameter

3. Numerical Experiments t ¡ t 0 ∞ ∞ + 1/Ø AEØ (a) Rectangular kernel ¡ t = AEØ [0,1/Ø] (t °∞) log t log ¡ t °log Ø log AEØ∞ slope º °(1 + ∞) (b)
Power law kernel on log-log scale for the rectangular one. We consider a non-symmetric block-matrix G to show that our method can e ectively uncover causality between the nodes, see Figure III.2. The matrix G has constant entries AE on the three blocks -AE = g i j = 1/6 for dimension 10 and AE = g i j = 1/10 for dimension 100 -, and zero outside. The two other parameters' values are the same for dimensions 10 and 100. The parameter ∞ is set to 1/2 on the three blocks as well, but we set three very di erent Ø 0 , Ø 1 and Ø 2 from one block to the other, with ratio

¡ t = AEØ∞(1 + Øt ) °(1+∞) t ¡ t 0 1/Ø AEØ (c) Exponential kernel ¡ t = AEØ exp(°Øt )
Ø i +1 /Ø i = 10
and Ø 0 = 0.1. The number of events is roughly equal to 10 5 on average over the nodes. We ran the algorithm on three simulated datasets: a 10-dimensional process with rectangular kernels named Rect10, a 10-dimensional process with power law kernels named PLaw10 and a 100-dimensional process with exponential kernels named Exp100.

MemeTracker dataset. We use events of the most active sites from the MemeTracker dataset2 . This dataset contains the publication times of articles in many websites/blogs from August 2008 to April 2009, and hyperlinks between posts. We extract the top 100 media sites with the largest number of documents, with about 7 million of events. We use the links to trace the Øow of information and establish an estimated ground truth for the matrix G. Indeed, when an hyperlink j appears in a post in website i , the link j can be regarded as a direct ancestor of the event. Then, Eq. (2) shows g i j can be estimated by

N i √ j T /N j T = #{links j ! i }/N j T .
Order book dynamics. We apply our method to AEnancial data, in order to understand the self and cross-inØuencing dynamics of all event types in an order book. An order book is a list of buy and sell orders for a speciAEc AEnancial instrument, the list being updated in real-time throughout the day. This model has AErst been introduced in [BJM16], and models the order book via the following 8-dimensional point process:

N t = (P (a) t , P (b) t , T (a) t , T (b) t , L (a) t , L (b) t ,C (a) t ,C (b) t ),
where P (a) (resp. P (b) ) counts the number of upward (resp. downward) price moves, T (a) (resp.

T (b)
) counts the number of market orders at the ask 3 (resp. at the bid) that do not move the price, L (a) (resp. L (b) ) counts the number of limit orders at the ask 4 (resp. at the bid) that do not move the price, and C (a) (resp. C (b) ) counts the number of cancel orders at the ask5 (resp. at the bid) that do not move the price. The AEnancial data has been provided by QuantHouse EUROPE/ASIA, and consists of DAX future contracts between 01/01/2014 and 03/01/2014.

Baselines. We compare NPHC to state-of-the art baselines: the ODE-based algorithm (ODE) by [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], the Granger Causality-based algorithm (GC) by [START_REF] Xu | Learning granger causality for hawkes processes[END_REF], the ADM4 algorithm (ADM4) by [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], and the Wiener-Hopf-based algorithm (WH) by [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF].

Metrics. We evaluate the performance of the proposed methods using the computing time, the Relative Error

RelErr(A,B) = 1 d 2 X i , j |a i j °bi j | |a i j | {a i j 6 =0} + |b i j | {a i j =0}
and the Mean Kendall Rank Correlation

MRankCorr(A,B) = 1 d d X i =1 RankCorr([a i • ], [b i • ]),
where RankCorr(x, y) = 2 d (d °1) (N concordant (x, y) °Ndiscordant (x, y)) with N concordant (x, y) the number of pairs (i , j ) satisfying x i > x j and y i > y j or x i < x j and y i < y j and N discordant (x, y) the number of pairs (i , j ) for which the same condition is not satisAEed.

Note that RankCorr score is a value between °1 and 1, representing rank matching, but can take smaller values (in absolute value) if the entries of the vectors are not distinct. Discussion. We perform the ADM4 estimation, with exponential kernel, by giving the exact value Ø = Ø 0 of one block. Let us stress that this helps a lot this baseline, in comparison to NPHC where nothing is speciAEed on the shape of the kernel functions. We used M = 10 basis functions for both ODE and GC algorithms, and L = 50 quadrature points for WH. We did not run WH on the 100-dimensional datasets, for computing time reasons, because its complexity scales with d 4 . We ran multi-processed versions of the baseline methods on 56 cores, to decrease the computing time.

Our method consistently performs better than all baselines, on the three synthetic datasets, on MemeTracker and on the AEnancial dataset, both in terms of Kendall rank correlation and estimation error. Moreover, we observe that our algorithm is roughly 50 times faster than all the considered baselines. On Rect10, PLaw10 and Exp100 our method gives very impressive results, despite the fact that it does not uses any prior shape on the kernel functions, while for instance the ADM4 baseline do. On Figure III.2, we observe that the matrix b G estimated with ADM4 recovers well the block for which Ø = Ø 0 , i.e. the value we gave to the method, but does not perform well on the two other blocks, while the matrix b G estimated with NPHC approximately reaches the true value for each of the three blocks. On these simulated datasets, NPHC obtains a comparable or slightly better Kendall rank correlation, but improves a lot the relative error.

On MemeTracker, the baseline methods obtain a high relative error between 9% and 19% while our method achieves a relative error of 7% which is a strong improvement. Moreover, NPHC reaches a much better Kendall rank correlation, which proves that it leads to a much better recovery of the relative order of estimated inØuences than all the baselines. Indeed, it has been shown in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF] that kernels of MemeTracker data are not exponential, nor power law. This partly explains why our approach behaves better.

On the AEnancial data, the estimated kernel norm matrix obtained via NPHC, see Figure III.3, gave some interpretable results (see also [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]):

1. Any 2 £ 2 sub-matrix with same kind of inputs (i.e. Prices changes, Trades, Limits or Cancels) is symmetric. This shows empirically that ask and bid have symmetric roles.

2. The prices are mostly cross-excited, which means that a price increase is very likely to be followed by a price decrease, and conversely. This is consistent with the wavy prices we observe on AEnancial markets.

3. The market, limit and cancel orders are strongly self-excited. This can be explained by the persistence of order Øows, and by the splitting of meta-orders into sequences of 4. Technical details smaller orders. Moreover, we observe that orders impact the price without changing it. For example, the increase of cancel orders at the bid causes downward price moves.

Technical details

We show in this section how to obtain the equations stated above, the estimators of the integrated cumulants and the scaling coe cient ∑ that appears in the objective function. We then prove the theorem of the paper.

Proof of Equation (8)

We denote ∫(z) the matrix

∫ i j (z) = L z ≥ t ! E(d N i u d N j u+t ) dud t ° §i § j ¥ ,
where L z ( f ) is the Laplace transform of f , and

√ t = P n∏1 ¡ (?n) t
, where ¡ (?n) t refers to the n th auto-convolution of ¡ t . Then we use the characterization of second-order statistics, AErst formulated in [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF] and fully generalized in [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF],

∫(z) = (I d + L °z (™))L(I d + L z (™)) > ,
where L i j = § i ± i j with ± i j the Kronecker symbol. Since I d + L z (™) = (I d °Lz (©)) °1, taking z = 0 in the previous equation gives > , which gives us the result since the entry (i , j ) of the last equation gives

∫(0) = (I d °G) °1L(I d °G> ) °1, C = RLR
C i j = P m § m R i m R j m .

Proof of Equation (9)

We start from [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF], cf. Eqs. (48) to (51), and group some terms:

K i jk = X m § m R i m R j m R km + X m R i m R j m X n § n R kn L 0 (√ mn ) + X m R i m R km X n § n R j n L 0 (√ mn ) + X m R j m R km X n § n R i n L 0 (√ mn ).
Using the relations L

0 (√ mn ) = R mn °±mn and C i j = P m § m R i m R j m
, proves Equation (9).

Integrated cumulants estimators

For H > 0 let us denote ¢

H N i t = N i t +H °N i t °H .
Let us AErst remark that, if one restricts the integration domain to (°H , H ) in Eqs. (4) and (5), one gets by permuting integrals and expectations:

§ i d t = E(d N i t ) C i j d t = E ≥ d N i t (¢ H N j t °2H § j ) ¥ K i jk d t = E ≥ d N i t (¢ H N j t °2H § j )(¢ H N k t °2H § k ) ¥ °d t § i E ≥ (¢ H N j t °2H § j )(¢ H N k t °2H § k ) ¥ .
The estimators ( 11) and ( 12) are then naturally obtained by replacing the expectations by their empirical counterparts, notably

E(d N i t f (t )) d t ! 1 T X ø2Z i f (ø).
For the estimator (13), we shall also notice that

E((¢ H N j t °2H § j )(¢ H N k t °2H § k )) = ZZ [°H ,H ] (t ) [°H ,H ] (t 0 )C j k t °t 0 d td t 0 = Z (2H °|t |) + C j k t d t.
We estimate the last integral with the remark above.

Choice of the scaling coe cient ∑

Following the theory of GMM, we denote m(X , µ) a function of the data, where X is distributed with respect to a distribution P µ 0 , which satisAEes the moment conditions g (µ) = E[m(X , µ)] = 0 if and only if µ = µ 0 , the parameter µ 0 being the ground truth. For x 1 ,..., x N observed copies of X , we denote b

g i (µ) = m(x i , µ), the usual choice of weighting matrix is c W N (µ) = 1 N P N i =1 b g i (µ) b g i (µ) >
, and the objective to minimize is then

√ 1 N N X i =1 b g i (µ) ! °c W N (µ 1 ) ¢ °1 √ 1 N N X i =1 b g i (µ) ! , (15)
where µ 1 is a constant vector. Instead of computing the inverse weighting matrix, we rather use its projection on {AEI d : AE 2 R}. It can be shown that the projection choses AE as the mean eigenvalue of c W N (µ 1 ). We can easily compute the sum of its eigenvalues:

Tr( c W N (µ 1 )) = 1 N N X i =1
Tr( b

g i (µ 1 ) b g i (µ 1 ) > ) = 1 N N X i =1
Tr( b

g i (µ 1 ) > b g i (µ 1 )) = 1 N N X i =1 || b g i (µ 1 )|| 2 2 .

Technical details

In our case, b 2 for the second. We compute the previous terms with R 1 = 0. All together, the objective function to minimize is

g (R) = h vec[ c K c °K c (R)], vec[ b C °C (R)] i > 2 R 2d 2 .
1 k c K c k 2 2 kK c (R) °c K c k 2 2 + 1 k b C k 2 2 kC (R) °b C k 2 2 . (16)
Dividing this function by

≥ 1/k c K c k 2 2 + 1/k b C k 2 2 ¥ °1, and setting ∑ = k c K c k 2 2 /(k c K c k 2 2 + k b C k 2
2 ), we obtaind the loss function given in Equation (10).

Proof of the Theorem

The main di erence with the usual Generalized Method of Moments, see [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF], relies in the relaxation of the moment conditions, since we have E[ b g T (µ 0 )] = m T 6 = 0. We adapt the proof of consistency given in [START_REF] Newey | Large sample estimation and hypothesis testing[END_REF].

We can relate the integral of the Hawkes process's kernels to the integrals of the cumulant densities, from [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF]. Our cumulant matching method would fall into the usual GMM framework if we could estimate -without bias -the integral of the covariance on R, and the integral of the skewness on R 2 . Unfortunately, we can't do that easily. We can however estimate without bias

R f T t C i j t d t and R f T t K i jk t d t with f T a compact supported function on [°H T , H T ]
that weakly converges to 1, with H T °! 1. In most cases we will take

f T t = [°H T ,H T ] (t ).
Denoting b C i j,(T ) the estimator of

R f T t C i j t d t, the term | E[ b C i j,(T ) ]°C i j | = | R f T t C i j t d t °C i j |
can be considered a proxy to the distance to the classical GMM. This distance has to go to zero to make the rest of GMM's proof work: the estimator b C i j,(T ) is then asymptotically unbiased towards C i j when T goes to inAEnity.

Notations

We observe the multivariate point process (N t ) on R + , with Z i the events of the i th component. We will often write covariance / skewness instead of integrated covariance / skewness. In the rest of the document, we use the following notations.

Hawkes kernels' integrals

G true = R © t d t = ( R ¡ i j t d t) i j = I d °(R true ) °1 Theoretical mean matrix L = diag( § 1 ,..., § d ) Theoretical covariance C = R true L(R true ) > Theoretical skewness K c = (K i i j ) i j = (R true ) Ø 2 C > + 2[R true Ø (C °Rtrue L)](R true ) > Filtering function f T ∏ 0 supp( f T ) Ω [°H T , H T ] F T = R f T s d s e f T t = f T °t Events sets Z i ,T,1 = Z i \ [H T , T + H T ] Z j ,T,2 = Z j \ [0, T + 2H T ]
Estimators of the mean b

§ i = N i T +H T °N i H T T e § j = N j T +2H T T +2H T Estimator of the covariance b C i j,(T ) = 1 T P ø2Z i ,T,1 °Pø 0 2Z j ,T,2 f ø 0 °ø °e § j F T ¢ Estimator of the skewness 6 b K i jk,(T ) = 1 T X ø2Z i ,T,1 √ X ø 0 2Z j ,T,2 f ø 0 °ø °e § j F T !√ X ø 00 2Z k,T,2 f ø 0 °ø °e § k F T ! °b § i T + 2H T X ø 0 2Z j ,T,2 √ X ø 00 2Z k,T,2 ( f T ? e f T ) ø 0 °ø00 °e § k (F T ) 2 ! GMM related notations µ = R and µ 0 = R true g 0 (µ) = vec ∑ C °RLR > K c °RØ 2 C > °2[R Ø (C °RL)]R > ∏ 2 R 2d 2 b g T (µ) = vec " b C (T ) °R b LR > c K c (T ) °RØ 2 ( b C (T ) ) > °2[R Ø ( b C (T ) °R b L)]R > # 2 R 2d 2 Q 0 (µ) = g 0 (µ) > W g 0 (µ) b Q T (µ) = b g T (µ) > c W T b g T (µ)

Consistency

First, let's remind a useful theorem for consistency in GMM from [START_REF] Newey | Large sample estimation and hypothesis testing[END_REF].

Theorem 2. If there is a function Q 0 (µ) such that (i ) Q 0 (µ) is uniquely maximized at µ 0 ; (i i ) £ is compact; (i i i ) Q 0 (µ) is continuous; (i v) b Q T (µ) converges uniformly in probability to Q 0 (µ), then b µ T = arg max b Q T (µ) P °! µ 0 .
We can now prove the consistency of our estimator.

µ = µ 0 , 2. µ 2 £, which is compact, 6 When f T t = [°H T ,H T ] (t ), we remind that ( f T ? e f T ) t = (2H T °|t |) +
. This leads to the estimator we showed in the article.

Technical details

the spectral radius of the kernel norm matrix satisAEes ||©||

§ < 1, 4. 8i , j , k 2 [d ], R f T u C i j u du ! R C i j u du and R f T u f T v K i jk u,v dud v ! R K i jk u,v dud v,
(µ) = [W 1/2 g 0 (µ)] > [W 1/2 g 0 (µ)] > 0 = Q 0 (µ 0 )
. Indeed, there exists a neighborhood N of µ 0 such that µ 2 N \{µ 0 } and g 0 (µ) 6 = 0 since g 0 (µ) is a polynom. Condition 2.1(ii) follows by

(i i ). Condition 2.1(iii) is satisAEed since Q 0 (µ) is a polynom. Con- dition 2.1(i v) is harder to prove. First, since b g T (µ) is a polynom of µ, we prove easily that E[sup µ2£ | b g T (µ)|] < 1.
Then, by £ compact, g 0 (µ) is bounded on £, and by the triangle and Cauchy-Schwarz inequalities, The estimator of L is unbiased so let's focus on the variance of b

| b Q T (µ) °Q0 (µ)| ∑ |( b g T (µ) °g0 (µ)) > c W T ( b g T (µ) °g0 (µ))| + |g 0 (µ) > ( c W T + c W > T )( b g T (µ) °g0 (µ))| + |g 0 (µ) > ( c W T °W )g 0 (µ)| ∑ k b g T (µ) °g0 (µ)k 2 k c W T k + 2kg 0 (µ)kk b g T (µ) °g0 (µ)kk c W T k + kg 0 (µ)k 2 k c W T °W k. To prove sup µ2£ | b Q T (µ) °Q0 (
L. E[( b § i ° §i ) 2 ] = E " µ 1 T Z T +H T H T (d N i t ° §i d t) ∂ 2 # = 1 T 2 Z T +H T H T Z T +H T H T E[(d N i t ° §i d t)(d N i t 0 ° §i d t 0 )] = 1 T 2 Z T +H T H T Z T +H T H T C i i t 0 °t d td t 0 ∑ 1 T 2 Z T +H T H T C i i d t = C i i T °! 0
By Markov inequality, we have just proved that k b

L °Lk P °! 0. Proof that k b C (T ) °C k P °! 0 First, let's remind that E( b C (T ) ) 6 = C . Indeed, E ≥ b C i j,(T ) ¥ = E µ 1 T Z T +H T H T d N i t Z T +2H T 0 d N j t 0 f t 0 °t °b § i e § j F T ∂ = E µ 1 T Z T +H T H T d N i t Z T +2H T °t °t d N j t +s f s ° §i § j F T ∂ + ≤ i j,T,H T F T = 1 T Z T +H T H T Z H T °HT f s E ≥ d N i t d N j t +s ° §i § j d s ¥ + ≤ i j,T,H T F T = Z f s C i j s d s + ≤ i j,T,H T F T Now, ≤ i j,T,H T = E ≥ § i § j °b § i e § j ¥ = °1 T 2 Z T +H T H T Z T +2H T 0 E ≥ d N i t d N j t 0 ° §i § j d td t 0 ¥ = °1 T 2 Z T +H T H T Z T +2H T 0 C i j t °t 0 d td t 0 = °1 T Z µ 1 + µ H T °|t | T ∂ °∂+ C i j t d t Since f satisAEes F T = o(T ), we have E( b C (T ) ) °! C . It remains now to prove that k b C (T ) °E( b C (T ) )k P °! 0.
Let's now focus on the variance of b

C i j,(T ) : V( b C i j,(T ) ) = E °( b C i j,(T ) ) 2 ¢ °E( b C i j,(T ) ) 2 . Now, E ≥ ( b C i j,(T ) ) 2 ¥ = E √ 1 T 2 X (ø,¥,ø 0 ,¥ 0 )2(Z i ,T,1 ) 2 £(Z j ,T,2 ) 2 ( f ø 0 °ø °F T /(T + 2H T ))( f ¥ 0 °¥ °F T /(T + 2H T )) ! = E µ 1 T 2 Z t ,s2[H T ,T +H T ] Z t 0 ,s 0 d N i t d N j t 0 d N i s d N j s 0 ( f t 0 °t °F T /(T + 2H T ))( f s 0 °s °F T /(T + 2H T )) ∂ = 1 T 2 Z t ,s2[H T ,T +H T ] Z t 0 ,s 0 2[0,T +2H T ] E ≥ d N i t d N j t 0 d N i s d N j s 0 ¥ • ( f t 0 °t °F T /(T + 2H T ))( f s 0 °s °F T /(T + 2H T ))
And,

E( b C i j,(T ) ) 2 = 1 T 2 Z t ,s2[H T ,T +H T ] Z t 0 ,s 0 2[0,T +2H T ] E ≥ d N i t d N j t 0 ¥ E ≥ d N i s d N j s 0 ¥ • ( f t 0 °t °F T /(T + 2H T ))( f s 0 °s °F T /(T + 2H T ))
Then, the variance involves the integration towards the di erence of moments µ r,s,t ,u °µr,s µ t ,u . Let's write it as a sum of cumulants, since cumulants density are integrable. µ r,s,t ,u °µr,s µ t ,u = ∑ r,s,t ,u + ∑ r,s,t ∑ u

[4] + ∑ r,s ∑ t ,u [3] + ∑ r,s ∑ t ∑ u [6] + ∑ r ∑ s ∑ t ∑ u °(∑ r,s + ∑ r ∑ s )(∑ t ,u + ∑ t ∑ u ) = ∑ r,s,t ,u + ∑ r,s,t ∑ u + ∑ u,r,s ∑ t + ∑ t ,u,r ∑ s + ∑ s,t ,u ∑ r + ∑ r,t ∑ s,u + ∑ r,u ∑ s,t + ∑ r,t ∑ s ∑ u + ∑ r,u ∑ s ∑ t + ∑ s,t ∑ r ∑ u + ∑ s,t ∑ r ∑ u
In the rest of the proof, we denote

a t = t 2[H T ,T +H T ] , b t = t 2[0,T +2H T ] , c t = t 2[°H T ,H T ] , g t = f t °1 T +2H T F T
Before starting the integration of each term, let's remark that:

1. ™ t = P n∏1 © (?n) t ∏ 0 since © t ∏ 0.
f t 0 °t d td t 0 = T F T b) R a t b t 0 g t 0 °t d td t 0 = 0 c) R a t b t 0 |g t 0 °t |d td t 0 ∑ 2T F T 4. 8t 2 R, a t (b ? e g ) t = 0
, where e g s = g °s .

Fourth cumulant

We want here to compute

R ∑ i , j ,i , j t ,
t 0 °t g s 0 °s | ∑ (|| f || 1 (1 + 2H T /T )) 2 ∑ 4|| f || 2 1 . Ø Ø Ø 1 T 2 Z ∑ i , j ,i , j t ,t 0 ,s,s 0 a t b t 0 a s b s 0 g t 0 °t g s 0 °s d td t 0 d sd s 0 Ø Ø Ø ∑ µ 2|| f || 1 T ∂ 2 Z d t a t Z d t 0 b t 0 Z d sa s Z d s 0 b s 0 M i ji j t 0 °t ,s°t ,s 0 °t ∑ µ 2|| f || 1 T ∂ 2 Z d t a t Z d t 0 b t 0 Z d sa s Z d w M i ji j t 0 °t ,s°t ,w ∑ µ 2|| f || 1 T ∂ 2 Z d t a t Z M i ji j u,v,w dud vd w ∑ 4|| f || 2 1 T M i ji j °! T !1 0 Third £ First
We have four terms, but only two di erent forms since the roles of (s, s 0

) and (t , t 0 ) are symmetric. First form

Z ∑ i , j ,i t ,t 0 ,s § j G t d t = § j T 2 Z ∑ i , j ,i t ,t 0 ,s a t b t 0 a s b s 0 g t 0 °t g s 0 °s d td t 0 d sd s 0 = § j T 2 Z ∑ i , j ,i t ,t 0 ,s a t b t 0 a s (b ? e g ) s g t 0 °t d td t 0 d s = 0 since a s (b ? e g ) s = 0 Second form Ø Ø Ø Z ∑ i , j , j t ,t 0 ,s 0 § i G t d t Ø Ø Ø = Ø Ø Ø § i T 2 Z ∑ i , j , j t ,t 0 ,s 0 a t b t 0 a s b s 0 g t 0 °t g s 0 °s d td t 0 d sd s 0 Ø Ø Ø = Ø Ø Ø § i T 2 Z ∑ i , j , j t ,t 0 ,s 0 a t b t 0 g t 0 °t b s 0 (a ? g ) s 0 d td t 0 d s 0 Ø Ø Ø ∑ § i T 2 2|| f || 1 Z d s 0 b s 0 (a ? |g |) s 0 Z d t a t Z d t 0 b t 0 K i j j t 0 °s0 ,t °s0 ∑ 4|| f || 1 K i j j § i F T T °! T !1 0 Second £ Second First form Ø Ø Ø Z ∑ i ,i t ,s ∑ j , j t 0 ,s 0 G t d t Ø Ø Ø ∑ 2|| f || 1 T 2 Z C i i t °sC j j t 0 °s0 a t b t 0 |g t 0 °t |a s b s 0 d td t 0 d sd s 0 ∑ 2|| f || 1 T 2 C i i C j j Z a t b t 0 |g t 0 °t |d td t 0 ∑ 4|| f || 1 C i i C j j F T T °! T !1 0 Second form Ø Ø Ø Z ∑ i , j t ,s 0 ∑ i , j t 0 ,s G t d t Ø Ø Ø ∑ 4|| f || 1 (C i j ) 2 F T T °! T !1 0 Second £ First £ First First form Z ∑ i , j t ,t 0 § i § j G t d t = § i § j T 2 Z ∑ i , j t ,t 0 a t b t 0 g t 0 °t d td t 0 Z a s b s 0 g s 0 °s d sd s 0 = 0 Second form Z ∑ i ,i t ,s § j § j G t d t = µ § j T ∂2 Z ∑ i ,i t ,s a t b t 0 g t 0 °t a s (b ? e g ) s d td t 0 d s = 0
We have just proved that V( b

C (T ) ) P °! 0. By Markov inequality, it ensures us that k b C (T ) °E( b C (T ) )k P °! 0, and AEnally that k b C (T ) °C k P °! 0. Á Proof that k c K c (T ) °K c k P °! 0
The scheme of the proof is similar to the previous one. The upper bounds of the integrals involve the same kind of terms, plus the new term (F T ) 2 /T that goes to zero thanks to the assumption 5 of the theorem.

Conclusion

In this paper, we introduce a simple nonparametric method (the NPHC algorithm) that leads to a fast and robust estimation of the matrix G of the kernel integrals of a Multivariate Hawkes process that encodes Granger causality between nodes. This method relies on the matching of the integrated order 2 and order 3 empirical cumulants, which represent the simplest set of global observables containing su cient information to recover the matrix G. Since this matrix fully accounts for the self-and cross-inØuences of the process nodes (that can represent agents or users in applications), our approach can naturally be used to quantify the degree of endogeneity of a system and to uncover the causality structure of a network.

By performing numerical experiments involving very di erent kernel shapes, we show that the baselines, involving either parametric or non-parametric approaches are very sensible to model misspeciAEcation, do not lead to accurate estimation, and are numerically expensive, while NPHC provides fast, robust and reliable results. This is conAErmed on the MemeTracker database, where we show that NPHC outperforms classical approaches based on EM algorithms or the Wiener-Hopf equations. Finally, the NPHC algorithm provided very satisfying results on AEnancial data, that are consistent with well-known stylized facts in AEnance.

Introduction

The previous approach based on the Generalized Method of Moments need the AErst three cumulants to obtain enough information from the data to recover the d 2 entries of G. Indeed, we want to recover d 2 independent coe cients -the entries of G -and the AErst two integrated cumulants give d + d (d + 1)/2 independent terms since the integrated covariance C is a symmetric matrix. Assuming the matrix G has a certain structure, we can get rid of the third order cumulant and design another estimation procedure using only the AErst two integrated cumulants. The advantage of such approach lies in the convexity of the related optimization problem, on the contrary to the minimization of L T from the previous chapter. The matrix we want to estimate minimize a simple criterion f convex, typically a norm, while being consistent with the AErst two empirical integrated cumulants.

Problem setting

We start from the relation between the integrated covariance C and the matrix R introduced in the previous chapter, from [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] and many other references:

C = RLR > .
Our purpose is still to approximate G = I °R°1 from the information encoded in the integrated cumulants. The previous equation in R admits a set of roots of that is isomorphocic to orthogonal group O n (R), and then:

G = I °L1/2 MC °1/2 with M 2 O n (R)
i.e.

L °1/2 (I °G)C 1/2 2 O n (R).
(1)

The previous expression only comes from the relation on the covariance. However, two classic assuptions on the Hawkes kernel norm matrix are not yet encoded. The AErst one concerns the positivity of the kernels, and then the positivity of their integrals:

g i j ∏ 0 for i , j 2 [d ].
Some variants of Hawkes processes allow the possibility of modeling inhibition through negative valued kernels [START_REF] Pernice | How structure determines correlations in neuronal networks[END_REF], with nonlinear Hawkes processes for instance [START_REF] Brémaud | Stability of nonlinear hawkes processes[END_REF], but the closed formulas of the cumulants [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] no longer stand with those variants. The other well-known assumption is linked to the stationarity of the process. The counting process N t has asymptotically stationary increments if the spectral norm of the kernel norm matrix is smaller than one: ||G|| < 1.

We AEnally encode the structure of the matrix G via the minimization of a criterion f subject to some constraints. For the problem to be easy to solve, the criterion f will be a convex function whose proximal operator is explicit. The AEt to the data encoded in Equation (1), and the two assumptions above will be regarded as constraints of our optimization problem.

All together, we formulate our problem as the following constrained optimization problem Constrained optimization problem:

min G f (G) s.t. L °1/2 (I °G)C 1/2 2 O n (R) ||G|| < 1 g i j ∏ 0
The problem above involves easy constraints on two di erent matrices: G and M = L °1/2 (I °G)C 1/2 . Our AErst idea is to relax the previous problem to turn it into a convex optimization problem.

The objective f is convex, and the constraints ||G|| < 1 and g i j ∏ 0 correspond to convex sets. The constraint that involves the orthogonal group is trickier and is not classic. We prove in Section 6 that the convex hull of the orthogonal group O n (R) is the closed unit ball w.r.t. the `2 norm. In the rest of the chapter, we denote B (resp. B

2 ) the open (resp. closed) unit ball w.r.t the spectral norm (resp. the `2 norm).

Instead of the previous problem, we split the variables G and M, meaning that we focus on the minimization problem both on G and M. Such minimization problem on two variables

x and z linked via an equation of the form Ax + B z = c can be e ciently solved with the Alternating Direction Method of Multipliers algorithm [START_REF] Glowinski | Sur l'approximation, par éléments AEnis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de dirichlet non linéaires[END_REF][START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via AEnite element approximation[END_REF] detailed in Section 3.

ADMM

The minimization problem we AEnally aim at solving writes:

min G,M f (G) + B 2 (M) + B (G) + R d £d + (G) (2) s.t. L °1/2 G + M C °1/2 = L °1/2
, On the contrary to the optimization problem of the previous chapter, the problem just stated is convex. We test this procedure on numerical simulations of various Hawkes kernels and real order book data, and we show how the criterion f impact the matrices we retrieve.

ADMM

The ADMM algorithm

The Alternating Direction Methods of Multipliers (ADMM) is a widely-used minimization method to solve constrained problems of the form

min x,z f (x) + g (z) (3) s.t. Ax + B z = c.
The objective function is separable in (x, z) with g and h two convex functions. The constraint involves two matrices A and B , and a constant vector c. The algorithm ADMM was originally introduced in [GM76] and [START_REF] Glowinski | Sur l'approximation, par éléments AEnis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de dirichlet non linéaires[END_REF], and focuses on the augmented Lagrangian [START_REF] Hestenes | Multiplier and gradient methods[END_REF][START_REF] Powell | A method for non-linear constraints in minimization problems[END_REF] associated to problem (3), that is:

L Ω (x, z, y) := g (x) + h(z) + y > (Ax + B z °c) + Ω 2 ||Ax + B z °c|| 2 2 , (4)
with Ω > 0 and solves the problem

min x,z max y L Ω (x, z, y) (5)
instead of the initial one. The method of multipliers [START_REF] Hestenes | Multiplier and gradient methods[END_REF][START_REF] Powell | A method for non-linear constraints in minimization problems[END_REF] (analysis in [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF]) applied to this problem would alternate an exact minimization step on the primal variable (x, z) and a gradient ascent step on the dual variable y. Instead of the exact minimization step on the couple (x, z), we do one pass of a Gauss-Seidel method [START_REF] Golub | Matrix computations[END_REF] and split the joint minimization into two partial minimization steps: one over x with z AExed, the other over z with x AExed. These two minimization steps can be done simultaneously, from the same initial points, or in the case of ADMM, one after the other, with an update between. Namely, ADMM algorithm iterates the following update steps:

x t +1 = argmin x L Ω (x, z t , y t ), z t +1 = argmin z L Ω (x t +1
, z, y t ),

y t +1 = y t + Ω(Ax t +1 + B z t +1 °c).

Convergence results

The convergence results of ADMM hold under the following two assumptions:

• The functions f and g are convex, proper1 and closed2 .

• The (unaugmented) Lagrangian L 0 has a saddle point i.e. there exist (x § , z § , y §

) for which

L 0 (x § , z § , y) ∑ L 0 (x § , z § , y § ) ∑ L 0 (x, z, y § )
for all x, y, z.

Under these two assumptions, the ADMM iterates satisfy the following convergences (a proof is given in [BPC + 11]):

• Residual convergence: r t = Ax t + B z t °c ! 0 as t ! 1 i.e. the iterates approaches feasibility.

• Objective convergence:

f (x t ) + g (z t ) ! min x,z { f (x) + g (z)} as t ! 1 i.e.
the objective function approaches its optimal value.

• Dual variable convergence: y t ! y § as t ! 1, where y § is a dual optimal point.

Examples

The ADMM method is quite general and plenty of optimization problems can be solved with it. We show here two usual tricks to turn an optimization problem into a relevant ADMM form. The AErst is to introduce indicator functions and concerns for instance optimization problem constrained on a set C :

min x f (x) s.t. x 2 C .
This problem can be equivalenty written:

min x,z f (x) + g (z)
s.t. x °z = 0 with g to be indicator of C i.e. to equal zero on C and 1 outside.

The other trick is to introduce a variable z being equal to a linear transformation of x. We consider the problem called total variation denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]:

min x ||x °b|| 2 2 + ∏ d °1 X i =1 |x i +1 °xi |.

Numerical results

Denoting F = (F i j ) with F i j = 1 j =i +1 °1j=i , the previous problem can be written as:

min x,z ||x °b|| 2 2 + ∏||z|| 1 s.t. F x °z = 0.
Such problem can be e ciently solved using the ADMM algorithm, since each update step of the algorithm has a closed form using proximal operators of the `2 and `1 norms.

Numerical results

In the previous sections, we only assumed f was a convex criterion whose proximal operator can be easily computed. Now, we exhibit three di erent choices for f and present the results obtained with these choices for both simulated and real-world dataset. The criteria we consider are the `1-norm

f = || • || 1 , Problem II when f = || • || 2
2 and Problem III for the case f = || • || § . We solve those minimization problem using the ADMM algorithm whose update steps are written above. The explicit update steps are provided in Section 6.

Simulated data

We simulated multivariate Hawkes point processes with the procedure already explained in the previous chapter, and implemented in the open-source library tick. As previously, we simulated three datasets generated from three di erent Hawkes kernels: the exponential kernel, the power-law kernel and the rectangular one. The mean vector and the integrated covariance matrix are computed using estimators provided in the previous chapter. We then used ADMM algorithm to solve the problems I , I I and I I I , and observed the same patterns for the three kernels. To ease the reading, we only show the results for the exponential kernel in dimension 100. The results from the Figure IV.1 and the Table IV.1 are consistent and shows that the solution to Problem I is the closest to the ground-truth matrix G. Moreover, according the Figure IV.1 one observes that the solutions to the two other problems are symmetric matrices, while the ground-truth matrix is not.

Order book data

The numerical experiments on simulated data incites us to focus on Problem I if the matrix G we want to uncover is not symmetric. Such non-symmetric relationships are for instance highlighted in [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF], where the authors studied the interplay between orders of di erent sizes. Indeed, a large trade is more likely to be followed by smaller trades than the opposite.

We use the same data as the authors of [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF] i.e. high-frequency order book data of futures traded at EUREX, that are also used in the numerical part of the previous chapter, see this section for details about the dataset. Here, we use the trades' timestamps of Bund futures.

We consider here unsigned trades i.e. we do not distinguish between buyer initiated trades and seller initiated ones. The di erent dimensions of the multivariate point process correspond to di erent intervals of volumes: each transaction falls into only one component. We denote N a t the number of transactions whose volume equals a that occured before t , N a:b t the number of transactions whose volume is between a and b (included), and N a: t the number of transactions whose volume is greater or equal than a. We then consider the following multivariate point processes, and solve the Problem I for these timestamped events: t . This solution is consistent with estimtates in lower dimension.

A t = (N 1 t ,
The solutions we found share the same patterns. We observe that self-excitation is preponderant, followed in importance by the excitation from large volumes. The excitation from large volumes is however lighter when we increase the dimension, this may be a consequence of the `1 norm minimization which aims at AEnding sparse solutions. Our observations are consistent with the results obtained in [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF], see this reference for AEnancial interpretations of the results. Note that kernel norm matrix of 21-dimensional model who have been very long to estimate with the Wieher-Hopf based method used in [RBL17], while our method has way lower complexity (comparable to the NPHC's one, see the previous chapter for a full comparison).

The approach consisting of minimizing a criterion instead of extracting the information from the third integrated cumulant seems promising. The convexity of the optimization problem is a real advantage compared to the non-convex problem one has to solve to estimate the Hawkes kernel norm matrix using NPHC, see Chapter III. The method developed in this chapter however lacks theory, especially for the choice of the criterion to minimize. As shown in the numerical part, the solution to the Problem I seems to provide better solutions, compared to the two other problems. Such statement could be explained by the similarity between Problem I and the `1 minimization of the compressing sensing problem

min x ||x|| 1 s.t. Ax = b.
One can indeed prove the exact recovery of the vector x under some assumptions [START_REF] Donoho | Compressed sensing[END_REF].

Technical details

Convex hull of the orthogonal group

The convex hull of the orthogonal group is the unit ball for the `2 norm. This is a nice exercise that can be solved using simple tools of linear algebra. A proof can be found in [START_REF] Giorgi | Mathematics of optimization: smooth and nonsmooth case[END_REF] for instance.

Updates of ADMM steps 6.2.1 Notations

We AErst denote the functions used in 2: f

1 (X ) = f (X ), f 2 (X ) = R d £d + (X ), f 3 (X ) = B (X ) and f 4 (X ) = B 2 (X ).
We also denote A = C = L °1/2 and B = C °1/2 . After splitting to the right number of variables (so that the update steps of ADMM algorithm for problem 2 write with closed formula), the problem 2 becomes:

min X 1 ,X 2 ,X 3 ,X 4 ,Y 1 ,Y 2 f 1 (X 1 ) + f 2 (X 2 ) + f 3 (X 3 ) + f 4 (X 4 ) s.t. Y 1 + Y 2 = C X 1 °X2 = 0 X 3 °X2 = 0 A °1Y 1 °X1 = 0 Y 2 B °1 °X4 = 0 90
6. Technical details

Update steps

Now, the update steps of ADMM algorithm (using the scaled dual form, see [BPC + 11]) write:

X t +1 1 = argmin X 1 f 1 (X 1 ) + (Ω/2)||X 1 °X t 2 +U t 2 || 2 F + (Ω/2)||A °1Y t 1 °X1 +U t 4 || 2 F X t +1 2 = argmin X 2 f 2 (X 2 ) + (Ω/2)||X t +1 1 °X2 +U t 2 || 2 F + (Ω/2)||X t 3 °X2 +U t 3 || 2 F X t +1 3 = argmin X 3 f 3 (X 3 ) + (Ω/2)||X 3 °X t +1 2 +U t 3 || 2 F X t +1 4 = argmin X 4 f 4 (X 4 ) + (Ω/2)||Y t 2 B °1 °X4 +U t 5 || 2 F Y t +1 1 = argmin Y 1 ||Y 1 + Y t 2 °C +U t 1 || 2 F + ||A °1Y 1 °X t +1 1 +U t 4 || 2 F Y t +1 2 = argmin Y 2 ||Y t +1 1 + Y 2 °C +U t 1 || 2 F + ||Y 2 B °1 °X t +1 4 +U t 5 || 2 F U t +1 1 = U t 1 + (Y t +1 1 + Y t +1 2 °C ) U t +1 2 = U t 2 + (X t +1 1 °X t +1 2 ) U t +1 3 = U t 3 + (X t +1 3 °X t +1 2 ) U t +1 4 = U t 4 + (A °1Y t +1 1 °X t +1 1 ) U t +1 5 = U t 5 + (Y t +1 2 B °1 °X t +1 4 )

Proximal operators

The previous update steps can be written using proximal operators of the functions f

1 , f 2 , f 3 and f 4 .
Proximal operator of f 2 This one is straightforward:

prox f 2 (X ) = (X ) + = (max(x i j , 0)) i j )
Proximal operator of f 3 Using techniques given in [START_REF] Boyd | Convex optimization[END_REF], one easily shows that First, it is easier to compute the proximal operator of f AE

3 = ae 1 (•)∑AE for AE < 1. Let's SVD some matrix X 2 R d £d : X = U SV > where U ,V 2 O d (R) and S = diag(ae 1 ,...,ae d ) is diagonal. Using techniques given in [BV04], one easily shows that prox f AE 3 (X ) = d X i =1 (ae i °(ae i °AE) + ) u i v > i Proximal operator of f 4
This one is a well-known projection too:

prox f 4 (X ) = X 1 {||X || 2 ∑1} + X ||X || 2 1 {||X || 2 >1} .

Final algorithm

X t +1 1 = prox f 1 /(2Ω) °(X t 2 °U t 2 + A °1Y t 1 +U t 4 )/2 ¢ X t +1 2 = (1/2) °X t +1 1 +U t 2 + X t 3 +U t 3 ¢ + X t +1 3 = prox f AE 3 °X t +1 2 °U t 3 ¢ X t +1 4 = prox f 4 °Y t 2 B °1 +U t 5 ¢ Y t +1 1 = (I d + A °2) °1( A °1(X t +1 1 °U t 4 ) °Y t 2 +C °U t 1 ) Y t +1 2 = ((X t +1 4 °U t 5 )B °1 °Y t +1 1 +C °U t 1 )(I d + B °2) °1 U t +1 1 = U t 1 + (Y t +1 1 + Y t +1 2 °C ) U t +1 2 = U t 2 + (X t +1 1 °X t +1 2 ) U t +1 3 = U t 3 + (X t +1 3 °X t +1 2 ) U t +1 4 = U t 4 + (A °1Y t +1 1 °X t +1 1 ) U t +1 5 = U t 5 + (Y t +1 2 B °1 °X t +1 4 )

Introduction

With the large number of empirical studies devoted to high frequency AEnance, relying on datasets of increasing size and quality, many progresses have been made during the last decade in the modelling and understanding the microstructure of AEnancial markets. Within this context, as evidenced by this special issue, Hawkes processes have become a very popular class of models. The main reason is that they allow one to account for the mutual inØuence of various types of events in a simple and parsimonious way through a conditional intensity vector. Hawkes processes have been involved in many di erent problems of high frequency AEnance ranging from the simple description of the temporal occurrence of market orders or price changes ([Bow07, HB14, FS12]), to the complex modelling of the arrival rates of various kinds of events in a full order book model ([Lar07, Tok11, JA13]). We refer to [START_REF] Bacry | Hawkes processes in AEnance[END_REF] for a recent review. A multivariate Hawkes model of dimension d is characterized by a d£d matrix of kernels, whose elements ¡ i j (t ) account for the inØuence, after a lag t , of events of type j on the arrival rate of events of type i . The challenging issue of the statistical estimation of the shape of these excitation kernels has been addressed by many authors and various solutions have been proposed whose performances (accuracy and computational complexity) strongly depend on the empirical situation one considers. Indeed, if non-parametric methods like e.g. the EM method ([LM11]), the Wiener-Hopf method ([BM14a, BM16, BJM16]) or the contrast function method ([RBRGTM14]) can be applied in low dimensional situations with a large number of events, one has to consider parametric penalized alternatives (like e.g., in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF][START_REF] Yang | Mixture of mutually exciting processes for viral di usion[END_REF]) when one has to handle a system of very large dimension with a relative low number of observed events (as, e.g., when studying events associated with the node activities of some social networks).

As far as (ultra) high frequency AEnance is concerned, the overall number of events can be very large. These events occur in a very correlated manner (with long-range correlations) and the system dimensionality can vary from low to moderately high. In a series of recent papers, Bacry et al. have shown that the non parametric Wiener-Hopf method provides reliable estimations in order to describe, within a multivariate Hawkes model, various aspects of level-I order book Øuctuations: the coupled dynamics of mid-price changes, market and limit order arrivals ([BM14a, BJM16]), the impact of market orders ( [START_REF] Bacry | Market impacts and the life cycle of investors orders[END_REF]) or the interplay between book orders of di erent sizes ([RBL17]). However, if one wants to account for systems of larger dimensionality by considering for instance a wider class of event types or the book events associated with a basket (e.g. a couple) of assets, then the Wiener-Hopf method (or any other similar non-parametric method) may reach its limits as respect to both computational cost and estimation accuracy. On the other hand, a parametric approach can lead to strong bias in the estimated inØuences between components.

For this reason, in the present paper, we propose to estimate Hawkes models of order book data using the faster and simpler non-parametric approach introduced in [ABG + 17]. This method focuses only on the global properties of the Hawkes process. More precisely, it aims at estimating directly the matrix of the kernel norms (also called the branching ratio matrix) without paying attention to the precise shape of these kernels. As recalled in the next section, this matrix does not bring all the information about the process dynamics, but is su cient to disentangle the complex interactions between various type of events and estimate the magnitude of their self-and cross-excitations. Moreover, it allows one to estimate the amplitude of Øuctuations of endogenous origin as compared to those of exogenous sources. The method we propose can be considered as the multivariate extension of the approach pioneered by [START_REF] Hardiman | Branching-ratio approximation for the self-exciting hawkes process[END_REF] that proposed to estimate the kernel norm of a one-dimensional Hawkes model directly from the integral of the empirical correlation function. Unfortunately their approach cannot be immediately extended to a multivariate framework because it does not bring a su cient number of constraints as compared to the number of unknown parameters. The method of [ABG + 17] circumvents this di culty by taking into account the AErst three integrated cumulant tensors of Hawkes process.

The paper is organized as follows: in Section 2 we provide the main deAEnitions and properties of multivariate Hawkes processes and we introduce the main notations we use all along the paper. The cumulant method of Achab et al. is described and illustrated in Section 3. In Section 4 we estimate the matrix of kernel of Hawkes models for level-I book 2. Hawkes processes: deAEnitions and properties events associated with 4 di erents very liquid assets, namely DAX, Euro-Stoxx, Bund and Bobl future contracts. We AErst consider the 8-dimensional model proposed in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] in order to compare our method to the former results obtained with a computationally more complex Wiener-Hopf method. We then show that the cumulant approach can easily be extended to a 12-dimensional model where all types of level-I book events are considered. Within this model, we uncover all the relationships between these types of events and we study the daily amplitude variations of exogenous intensities. In Section 5 we investigate the correlation between two assets by considering the events of their order book within a 16-dimensional model. This allows us to discuss the inØuence of both their tick size and their degree of reactivity with respect to the impact of their book events on each other. Section 6 contains concluding remarks while some technical details are provided in Appendix.

Hawkes processes: deAEnitions and properties

In this section we provide the main deAEnitions and properties of multivariate Hawkes processes and set the notations we need all along the paper.

Multivariate Hawkes processes and the branching ratio matrix G

A multivariate Hawkes process of dimension d is a d -dimensional counting processes N t with a conditional intensity vector ∏ t that is a linear function of past events. More precisely,

∏ i t = µ i + d X j =1 Z t °1 ¡ i j (t °s) d N j s (1)
where µ i represents the baseline intensity while the kernel ¡ i j (t ) quantiAEes the excitation rate of an event of type j on the arrival rate of events of type i after a time lag t . In general it is assumed that each kernel is causal and positive, meaning that Hawkes processes can only account for mutual excitation e ects since the occurrence of some event can only increase the future arrival intensity of other events. In order to consider the possibility of inhibition e ects, one can allow kernels to take negative values. In that case, we have to consider expression (1) only when it provides a positive result while the conditional intensity is assumed to be zero otherwise. Rigorously speaking, such non-linear variant of Eq. (1) cannot be handled as simply as the original Hawkes process ([BM96]) but, as empirically shown in e.g. [START_REF] Reynaud-Bouret | Goodness-of-AEt tests and nonparametric adaptive estimation for spike train analysis[END_REF] or [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF], if the probability that ∏ i t < 0 is small enough, one can safely consider the model as linear so that all standard expressions provide accurate results. In the following we will suppose that we are in this case and we don't necessarily impose that the kernels ¡ i j (t ) are positive functions.

Let us deAEne the matrix G as the matrix whose coe cients are the integrals of the kernels

¡ i j
(t ) (that are supported by R + ):

G i j = Z +1 0 ¡ i j (t )d t . (2)
Let us remark that, as it can directly be seen from the cluster representation of Hawkes processes ([HO74]), G i j represents the mean total number of events of type i directly triggered by an event of type j . For that reason, in the literature, the matrix G is also referred to as the branching ratio matrix ([HB14]). Notice that since the kernels ¡ i j (t ) are not necessarily non negative functions, G i j does not in general correspond to the L 1 norm of ¡ i j . For the sake of simplicity, though this is not technically correct, we shall often refer to the matrix G as the "matrix of kernel norms" or more simply the "norm matrix".

If kGk stands for the largest eigenvalue of G, it is well known that a su cient condition for the intensity process ∏ t to be stationary is that kGk < 1. In the following we will always consider this condition satisAEed. One can then deAEne the matrix R as:

R = (I d °G) °1, (3)
where I d denotes the identity matrix of dimension d .

Let § denote the mean intensity vector:

§ = E(∏ t ) , (4)
so that the ratio µ i § i represents the fraction of events of type i that are of exogenous origin. One can easily prove that § and µ are related as:

§ = R µ (5)
If one deAEnes the matrix ™ as:

™ = GR = R °Id , (6)
then ™ i j represents the average number of events of type i triggered (directly or indirectly) by an exogenous event of type j . When one analyzes empirical data within the framework of Hawkes processes, the previous remarks allow one to quantify causal relationships between events in the sense of Granger, i.e., within a well deAEned mathematical model. In that respect, the coe cients of the matrices G or ™ can be read as (Granger-)causality relationships between various types of events and used as a tool to disentangle the complexity of the observed Øow of events occurring in some experimental situations [START_REF] Eichler | Graphical modeling for multivariate hawkes processes with nonparametric link functions[END_REF]). Let us emphasize that such causal implications are just a matter of interpretation of data within a speciAEc model (namely a Hawkes model) and should simply be considered as a convenient and parsimonious way to represent that data. They should not, in any way, be understood as a "physical" causality reØecting their "real nature".

Integrated Cumulants of Hawkes Process

The NPHC algorithm developed in [ABG + 17] and described in Sec. 3 below, enables the direct estimation of the matrix G from a single or several realizations of the process. It relies on the computation of low order cumulant functions whose expressions are recalled below.

Given 1 ∑ i , j , k ∑ d , the AErst three integrated cumulants of the Hawkes process can be, thanks to stationarity, deAEned as follows:

§ i d t = E(d N i t ) (7)
C i j d t = Z ø2R ≥ E(d N i t d N j t +ø ) °E(d N i t )E(d N j t +ø ) ¥ (8) K i jk d t = ZZ ø,ø 0 2R 2 ≥ E(d N i t d N j t +ø d N k t +ø 0 ) + 2E(d N i t )E(d N j t +ø )E(d N k t +ø 0 ) °E(d N i t d N j t +ø )E(d N k t +ø 0 ) °E(d N i t d N k t +ø 0 )E(d N j t +ø ) °E(d N j t +ø d N k t +ø 0 )E(d N i t ) ¥ , (9)
where Eq. ( 7) is the mean intensity of the Hawkes process, the second-order cumulant (8) refers to the integrated covariance density matrix and the third-order cumulant (9) measures the skewness of N t . Using the martingale representation ([BM16]) or the Poisson cluster process representation ([JHR15]), one can obtain an explicit relationship between these integrated cumulants and the matrix R (and therefore the matrix G thanks to Eq. ( 3)). Some straightforward computations (see [ABG + 17]) lead to the following identities:

§ i = d X m=1 R i m µ m (10) C i j = d X m=1 § m R i m R j m ( 11
) K i jk = d X m=1 (R i m R j m C km + R i m C j m R km +C i m R j m R km °2 § m R i m R j m R km ). (12)
3 The NPHC method

In this section we brieØy recall the main lines of the recent non parametric method proposed in [ABG + 17] that leads to a fast and robust direct estimation of the branching ratio matrix G without estimating the shape of the kernel functions. This method is based on the remark that, as shown in [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] and as it can be seen in Eqs. (10), ( 11) and ( 12), the integrated cumulants of a Hawkes process can be explicitly written as functions of R. The NPHC method is a moment method that consists in directly exploiting these equations to recover R and thus G.

Estimation of the integrated cumulants

Let us AErst introduce explicit formulas to estimate the three moment-based quantities listed in the previous section, namely, §, C and K . In what follows, we assume there exists H > 0 such that the truncation from (°1, +1) to [°H , H ] of the domain of integration of the quantities appearing in Eqs. ( 8) and (9) introduces only a small error. This amounts to neglecting tail e ects in the covariance density and in the skewness density, and it corresponds to a good approximation if (i ) each kernel ¡ i j (t ) is essentially supported by [0, H ] and (i i ) the spectral norm kGk is less than 1.

In this case, given a realization of a stationary Hawkes process {N t : t 2 [0, T ]}, as shown in [ABG + 17], we can write the estimators of the AErst three cumulants (7), ( 8) and (9) as

b § i = 1 T X ø2Z i 1 = N i T T (13) b C i j = 1 T X ø2Z i ≥ N j ø+H °N j ø°H °2H b § j ¥ (14)
b

K i jk = 1 T X ø2Z i ≥ N j ø+H °N j ø°H °2H b § j ¥ • ≥ N k ø+H °N k ø°H °2H b § k ¥ °b § i T X ø2Z j X ø 0 2Z k (2H °|ø 0 °ø|) + + 4H 2 b § i b § j b § k . (15)
In practice, the AEltering parameter H is selected by (i ) computing estimates of the covariance density at several points t 1 , (i i ) assessing the characteristic time ø c after which the covariance density is negligible, and (i i i ) setting a multiple of ø c for H , for instance H = 5ø c .

The NPHC algorithm

The covariance C only provides d (d + 1)/2 independent coe cients and is therefore not su cient to uniquely identify the d 2 coe cients of the matrix G. In order to set a su cient number of constraints, the NPHC approach relies on using all the covariance C along with a restricted number of the (d 3 + 3d 2 + 2d )/6 third-order independent cumulant components, namely the d 2 coe cients K c = {K i i j } 1∑i , j ∑d . Thus, we deAEne the estimator of R as b R 2

argmin R L (R), where

L (R) = (1 °∑)kK c (R) °c K c k 2 2 + ∑kC (R) °b C k 2 2 , (16)
where k • k 2 stands for the Frobenius norm, while c K c and b C are the respective estimators of C and K c as deAEned in Equations ( 14), (15) above. It is noteworthy that the above mean square error approach can be seen as a particular instance of Generalized Method of Moments (GMM), see [START_REF] Hall | Generalized method of moments[END_REF], [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF]. Though this framework allows to determine the optimal weighting matrix involved in the loss function, in practice this approach is unusable, as the associated complexity is too high. Indeed, since we have d 2 parameters, this matrix has d 4 coe cients and GMM calls for computing its inverse leading to a O(d 6 ) complexity. Thus, instead, we choose to use the loss function (16) in which, so as to be of the same order, the two terms are rescaled using ∑ = k c

K c k 2 2 /(k c K c k 2 2 +k b C k 2
2 ). We refer to Appendix 1 for an explanation of how ∑ is related to the weighting matrix. Finally the estimator of G is straightforwardly obtained as

b G = I d °b R °1,
from the inversion of Eq. (2). The authors of [ABG + 17] proved the consistency of the soobtained estimator b G, i.e. the convergence in probability to the true value, when the observation time T goes to inAEnity.

1 the pointwise covariance density at t can be estimated with Let us mention that, when applied to AEnancial time-series, the number of events is generally large as compared with d (i.e., n = max i |Z i | ¿ d ), thus the matrix inversion in the previous formula is not the bottleneck of the algorithm. Indeed, it has a complexity O(d 3 ) which is cheap as compared with the computation of the cumulants which is O(nd 2

1 hT P ø2Z i ≥ N j ø+t +h °N j ø+t °h b § j ¥ for a small h 102 3. The NPHC method t ¡ t 0 ∞ ∞ + 1/Ø AEØ (a) Rectangular kernel log t log ¡ t °log Ø log AEØ∞ slope º °(1 + ∞) (b)
). Thus, assuming the loss function ( 16) is minimized after N iter iterations, the overall complexity of the algorithm is O(nd 2 + N iter d 3

). The authors of [ABG + 17] compared the complexity of their algorithm with other state-of-the-art methods' ones, namely the ordinary di erential equations based (ODE) algorithm in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], the Sum of Gaussians based algorithm in [START_REF] Xu | Learning granger causality for hawkes processes[END_REF], the ADM4 algorithm in [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], and the Wiener-Hopf-based algorithm in [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF]. The complexity of NPHC is smaller, because the algorithm NPHC directly estimates the kernels' integrals while other methods go through the estimation of the kernel functions themselves.

Numerical experiments

As mentioned above, the NPHC algorithm is non parametric and provides an estimation of the integral of the kernels regardless of their shapes. In order to illustrate the stability of our method with respect to the shape of the kernels, we simulated two datasets with Ogata's Thinning algorithm introduced in [START_REF] Ogata | On lewis' simulation method for point processes[END_REF] using the open-source library tick2 . Each dataset corresponds to a di erent kernel shape (but with the same norm), a rectangular kernel and a power-law kernel, both represented in Figure V.1:

rectangular kernel: ¡(t ) = AEØ [0,1/Ø] (t °∞) (17) power law kernel: ¡(t ) = AEØ∞(1 + Øt ) °(1+∞) (18)
In both cases, AE corresponds the integral of the kernel, 1/Ø can be regarded as a characteristic time-scale, and ∞ corresponds to the scaling exponent for the power law kernel and a delay parameter for the rectangular one. We consider a non-symmetric 10-dimensional blockmatrix G with 3 non-zero blocks, and where the parameters AE = 1/6 and ∞ = 1/2 take the

T + T °L+ L °C + C °T a T b L a L b C a C b
DAX 11.9 11.9 21.8 21.9 10.1 10.1 11.6 11.7 80.0 79.5 97.3 96.1 ESXX 2.6 2.6 3.5 3.6 0.9 0.9 16.4 16.5 176.0 174.7 172.4 170.8 Bund 3.2 3.2 4.0 4.0 0.8 0.8 14.5 14.7 125.4 125.0 111.5 110.7 Bobl 1.1 1.1 1.5 1.5 0.5 0.5 6.1 6.1 86.5 86.8 81.6 81.4

Table V.1: Average number of events in thousands per type in a trading day (from open at 08:00 to closing at 22:00 Frankfurt time) for the four assets considered.

order and a trade is generated. It is therefore possible to obtain a list of the orders that were submitted complete with their time, type (limit, cancel or market order), volume and price. The timestamp precision is one microsecond and the timestamps are set directly by the exchange.

In this work we are interested in disentangling the interactions of di erent types of events occurring at the AErst level of the order book. To this end, we will distinguish the following event types:

• T + (T °) : upwards (downwards) mid price movement triggered by a market order;

• L + (L °) : upwards (downwards) mid price movement triggered by a limit order;

• C + (C °) : upwards (downwards) mid price movement triggered by a cancel order;

• T a (T b ) : market order at the ask (bid) that does not move the mid price;

• L a (L b ) : limit order at the ask (bid) that does not move the mid price;

• C a (C b ) : cancellation order at the ask (bid) that does not move the mid price.

Additionally, we introduce the symbols P + (P °) to denote an upwards (downwards) mid price movement irrespectively of its origin. In Table V.1 we report the average number of events per day (from 08:00 am to 10:00 pm) for each asset and each type. We remark that all four assets are extremely active securities with an average of more than 300.000 events per day.

One characteristic that strongly inØuences the order book dynamics at short time scales is the tick size to average spread ratio. When this ratio is close to one (resp. much smaller than one), the asset is said to be a "large tick asset" (resp. a "small tick asset") (see, e.g., [START_REF] Dayri | Large tick assets: implicit spread and optimal tick size[END_REF]). In our dataset, all assets are large-tick assets (the spread is equal to one tick in more than 95% of the times) except for the DAX future, which is a small-tick one. As evidenced by Table V.1, the price changes much less frequently on large tick assets. One can also remark that the quantity available at the best quotes tends do be proportionally much larger on large tick assets. These microstructural characteristics will be reØected by our analysis. method when the focus is solely on the kernel interaction matrix. Indeed, in order to estimate the kernel norm matrix with the Wiener-Hopf method, the full kernel functions have to be estimated AErst and then numerically integrated. The NPHC method thus represents a much faster alternative, as it does not require the estimation of d 2 functions but directly estimates their integrals. Besides the speed gain, the gain in complexity allows NPHC to scale much better when increasing the dimension, i.e., when using more detailed models.

A 12-dimensional mono-asset model

By estimating directly the norm of the kernels and not the whole kernel function, the NPHC method can be used to investigate systems of greater dimension. In this section we extend the model of Section 4.2 to 12 dimensions by separating the type of events that lead to a price move. The 12 even types we consider are thus

T + (T °), L + (L °), C + (C °), T a (T b ), L a (L b ), C a (C b ).
We then apply the NPHC algorithm to estimate the branching ratio matrix. When not otherwise speciAEed, we set H = 500s. To further assess the validity of our methodology and the impact of time-of-day e ects, we AErst estimate the model using di erent time slots within the trading day. In Section 4.3.2 we also check the robustness of our results as respect to the choice of the parameter H .

Kernel stability during the trading day

We ran our method for the DAX future on the 12-dimensional point process detailed above on di erent subintervals of the trading day. More precisely, we divided each trading day into 7 slots with edges at 08:00 am, 10:00 am, 12:00 am, 02:00 pm, 04:00 pm, 06:00 pm and 10:00 pm. We then estimated the 12-dimensional model described above on each slot separately, averaging over all 338 trading days available in our dataset. The results are remarkable in that the kernel norm matrix appears to be very stable during the trading day. (we checked that this is also true if we set H = 1s).

The NPHC method outputs the estimated matrix b R (and then b G) from which one can obtain an estimate of µ using the relation (5) that links R and the mean intensity §, namely

b µ = b R °1 b §.
In the right panel of Figure V.5 we plot the values of b µ as obtained using the above relation for the T a/b , L a/b and C a/b components. We consider the kernel norm matrix as constant in each two hours slot and we estimate the average intensity on 15 minutes non-overlapping windows. Moreover, for each type of events we show the average of the bid/ask components. For comparison, in the left panel of Figure V.5 we show the empirical intraday pattern obtained for each component. We remark that the values of µ obtained with our procedure vary during the day and roughly follow the intraday curve of the respective components. Let us notice that µ i / § i , the fraction of exogenous events, is of the order of a few percent. This is fully consistent with what was found in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] and means, within the Hawkes framework, that most of the observed order book dynamics is strongly endogenous. For the price moving components the values of § are of the order of 1s °1, while results for µ are more noisy, similarly to those of T a/b . This analysis conAErms the result formerly observed in [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF] that the kernels are stable during the day, and that time-of-day e ects are well captured by the baseline intensity, at least as long as we are mainly concerned with the high frequency dynamics on a very liquid asset as is the case here.

Analysis of the G matrix: Unveiling mutual interactions between book events

Having established that the estimated kernel matrix is stable with respect to time of the day e ects, we now examine more in-depth its structure. In Figure V.6 is represented the result of the estimation of the matrix G over the whole trading day for the DAX future. The branching ratio matrix on the left panel is estimated with H = 1s while the right panel corresponds to 4. Single-asset model H = 500s. Let us recall that both horizons are several orders of magnitude larger than the typical inter-event time.

Concerning the di erences between the two matrices, we note that certain inhibitory e ects that are visible for H = 1s are less intense or disappear when H = 500s is used. This most notably happens for the elements T + ! T + and T + ! T °and similarly for L + ! L +/°a nd C + ! C +/°, which suggests that when we look at longer scale correlation the self-exciting behavior (i.e. trades are followed by more trades) tends to prevail on the high frequency mean reverting e ect.

Apart from these di erences, we can make some observations that are valid in both cases. In particular, we note that two main interaction blocks stand out. The AErst is the upper left corner which concerns interactions between price-moving events, where two anti-diagonal bands are prominent. The second is the bottom right corner, which has a strong diagonal structure. The blocks involving interactions between price-moving and non-price moving events present instead much smaller values. In what follows, we AErst discuss more in depth the e ects of price movements on other events, then those of non-price-moving ones. We also remark that the spectral norm of the estimated matrices G is close to 1 while being inferior (e.g. 0.98 for the DAX with H = 500s). This is in line with what was found in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF], and the criticality of AEnancial markets highlighted in [START_REF] Hardiman | Critical reØexivity in AEnancial markets: a hawkes process analysis[END_REF].

Before entering into more details, let us remark that in both cases the expected symmetry up/down (+/°) and bid-ask (b/a) is well recovered in our results. Therefore, to make notation lighter and facilitate the exposition, we will comment only on one side. More precisely, when discussing the e ects of price moves we will refer only to the upwards ones (T + /L + /C + ) and when discussing e ects of liquidity changes we will focus on ask side events (T a /L a /C a ). E ect of price-moving events As we noted above, the most relevant interactions involving T + are the T + ! L + and T + ! L °ones, the mean reverting one (T + ! L °) being more intense. When a market order consumes the liquidity available at the best ask, two main scenarios can occur for the mid price to change again, either the consumed liquidity is replaced, reverting back the price (mean-reverting scenario, highly probable) or the price moves up again and a new best bid is created.

Market orders that move the mid price have also an inhibitory e ects at short time scales on subsequent price-moving trades (T + ! T + is negative for H = 1s). Indeed, once a market order consumes the liquidity available at the best quote, it is unlikely that the price will be moved in the same direction by other market orders as the price becomes more unfavorable. We also note a generally inhibitory e ect of T + on price-moving cancel orders which can be linked to a mechanical e ect, liquidity that has been consumed by the market order cannot be canceled anymore.

The same kind of dynamics is at play also in the interactions L + ! T + and L + ! T °with the roles inverted. Again, the mean reverting e ect L + ! T °appears to be much more probable. A strong mean-reverting e ect is found in the block L + ! C °. This is possibly the signature of high-frequency strategies whereby agents place limit orders in the spread and cancel them shortly thereafter.

Concerning C + events, the main feature lies in the block C + ! L °, where we notice the same anti-diagonal dominance found for the block L + ! C °. Again, we can suppose that when a limit order in the spread is removed it is often quickly replaced by market participants.

Finally, the e ect of price moving events on non-price moving ones can be summarized in two main e ects. The AErst is a trend-following/order splitting e ect by which e.g. trades at the ask are likely to be followed by more trades in the same direction (T + ! T a ) and similarly for limit (L + ! L b ) and cancel (C + ! C a ) orders. The second is the shift in liquidity triggered by a price change. A trade at the ask that moves upward the mid price triggers limit orders on the opposite side (T + ! L b ). This can be understood using a latent price argument ([RR10]), as it is well known that there are more limit orders far from the latent price. Right after the mid price goes up, the latent price is expected to be closer to the newly best ask price than to the best bid price, thus limit order Øow is expected to be higher at best bid than at best ask.

E ect of non-price-moving events

For all events T a , L a and C a the most visible feature is the strong self-exciting interaction. This has been conAErmed in several works ([BJM16], [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF]) and can be traced to order-splitting strategies and herding behaviors. Signatures of typical trading patterns can be seen also in the kernels L a ! C a , L a ! C b , where the positive value of the kernel arises form agents canceling and replacing their limit orders with or without switching sides.

We also note the positive e ects T a ! T + , L a ! T °and C a ! T + . All these e ects, as well as the analogous ones on C +/°/ L +/°, reØect the fact that changes in the imbalance have an inØuence on the probability of a subsequent price move. So when the queue at the best ask decrease an upward price move becomes more likely and vice-versa. These e ects are much more relevant on a small tick asset (DAX) than on a large tick asset (Bund) where, the size of the queues being larger, their inØuence is marginal.

We performed the same analysis on the Bund (see Figure V.7). The main di erences as compared to the DAX are that the e ects between events that move the price are much more intense while the e ects of events that do not move the price on those that do move the price (and vice-versa) are much less pronounced, indeed they are barely visible in Figure V.7. This can be basically seen as a simple consequence of the Bund future being large tick assets, while the DAX is a small tick one. Therefore, price movements on the former are much less frequent but when they happen their e ects are more marked.

Analysis of the ™ matrix: the AEngerprint of meta-orders

As discussed in Section 2, the elements of the matrix ™ quantiAEes the total e ect, direct and indirect, of an event of type j on events of type i . More precisely, thanks to the branching process structure, we can interpret √ i j as the mean number of events of type i generated by a single exogenous ancestor of type j . We plot the estimated matrices ™ for the DAX and . We note that an exogenous limit or cancel event generates a large number of limit and cancel events and, to a lesser extent, trade events. This can be read as the signature of meta-orders. Indeed, if an agent wants to sell a large number of contracts 6 , he will place a meta-order, i.e., he will optimize the overall cost by dividing this large order into several smaller orders. The overall optimization will result in many limit/cancel sell orders L a ,C a and, as less as possible, of sell market orders T b (the cost of a market order is on average higher than that of a limit order). The same description can be applied to understand why an exogenous sell market order T b generates mainly limit and cancel sell orders L a ,C a as well as other sell market orders T b .

Due to the much lower values of the exogenous intensities for price moving events, the left part of the ™ matrix is more noisy. Nevertheless, at least in the DAX case, we note also for the price moving components the prevalence of the L + ! L + and L + ! C °elements, which are the price-moving counterparts of the e ect described for L a .

Finally, we also remark that although we noted several inhibition e ects in the matrices G, the elements of ™ are non negative. This suggests that most inhibition e ects are short lived and the e ect of an event arrival is towards an increase of the overall intensity. This is in line with what was found in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] and [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF], where the inhibitions e ects were shown to be mostly concentrated around the typical market reaction time.

Within the branching ratio representation of Hawkes processes, µ j § i √ i j represents the fraction of events of type i that has a type j as primary ancestor. Along the same line, we can estimate the fraction of aggressive orders (i.e. all T ), as opposed to passive orders (L or 6 Let us recall that, in our discussion, we only address half of the matrix coe cients since the discussion on the other half can be obtained using the symmetries ask/bid, buy/sell, price up/price down. Following these lines, we only consider here the case of a selling meta-order. 112 5. Multi-asset model C ), that is ultimately generated by another aggressive order, as:

1 P i ={T +/°, T a/b } § i X j ={T +/°, T a/b } X i ={T +/°, T a/b } √ i j µ j . (19)
We AEnd that for both assets this fraction is about 10%, which means that the large majority of market orders have a "passive order" (L or C ) oldest ancestor. We compute the analogous fraction for passive orders and we AEnd that for both assets more than 96% of the passive orders (L or C ) have an oldest ancestor that is itself a passive (L or C ) order. This fact is in line with the idea that meta-orders would be at the origin of most of the trading activity within the order book.

Multi-asset model

Studying and quantifying the interactions and comovements within a basket of assets is an important topic in AEnance. Most of these studies focus on the return correlations properties in relationship with portfolio theory. At very high frequency, the discrete nature of price variations and the asynchronous occurrence of price change events make the correlation analysis trickier and, in order to avoid well known bias (like the Epps e ect) one has to use speciAEc techniques like the estimator proposed by [HY + 05]. Hawkes processes, being naturally deAEned in continuous time, can represent a complementary tool for the investigation of highfrequency cross-asset dynamics.

The idea of capturing the joint dynamic of multiple assets via Hawkes processes has only been considered in few recent papers. Let us mention the work proposed by [BCT + 15] which models the simultaneous cojumps of di erent assets using a one-dimensional Hawkes process, and a more recent work ( [START_REF] Fonseca | Correlation and lead-lag relationships in a hawkes microstructure model[END_REF]) which focuses on the correlation and lead-lag relationships between the price changes of two assets, in the spirit of [START_REF] Bacry | Modelling microstructure noise with mutually exciting point processes[END_REF].

In this section, we aim at unveiling a more precise structure of the high-frequency crossasset dynamics by pushing further the dimensionality of the model to include simultaneously events on two assets. We AErst consider the pair DAX-EURO STOXX and then the one Bobl-Bund. The pairs of assets considered here are tightly related, as they share exposure to the same risk factors and, in the case of DAX-EURO STOXX, also because the underlying indices actually share a signiAEcant part of their components. This is conAErmed also by Table V.2 where we report 5 minutes return correlations among the considered assets.

In this section we consider the same kind of events as in Section 4.2 and we have therefore a 16-dimensional model (2 £ 8) corresponding to 256 possible interactions. Let us point out that this is quite a large dimension value for a non parametric methodology.

The DAX -EURO STOXX model

In the following, we will denote the events of the DAX order book with the subscript D while we will use the subscript X for the events of EURO STOXX order book. The obtained branching ratio matrix is displayed in Figure V.9. We observe that the mono-asset submatrices (the two 8 £ 8 block matrices along the diagonal), which present the most relevant e ects, have the same structure as the ones which have already been commented on in detail in Section 4.2. Consequently, in this section, we shall focus our discussion on the non diagonal 8 £ 8 submatrices that correspond to the interactions between the two assets. These two submatrices are shown in Figure V.10. Note that colors have been rescaled to highlight their structure. To keep the notation lighter, we will comment only on e ect of upwards price moves and ask events as it was done in the previous section, since we AEnd the symmetries +/°and a/b to be well respected. The most striking feature emerging from Figure V.10 is the very intense relation between same-sign price movements on the two assets. Albeit present in both directions, the norms P + X ! P + D attain larger values. Another notable aspect is the di erent e ects of price moves and liquidity changes of one asset on events on the other asset. Price moves on the DAX have also an e ect on the Øow of limit orders on EURO STOXX (P + D ! L b X and P + D ! C a X ), whereas EURO STOXX price moves triggers mainly DAX price moves in the same direction (P + X ! P + D ). An important aspect for understanding this result is the di erent perceived tick sizes on the two assets.

In the following, whenever it is convenient, we shall place the discussion withing the framework of latent price models (e.g., [START_REF] Robert | A new approach for the dynamics of ultrahigh-frequency data: The model with uncertainty zones[END_REF]). Within this framework, the latent price refers to an underlying e cient price representing at any time some average opinion of market participants about the value of the asset. As noted in Section 4.1, the DAX future is a smalltick asset, while the EURO STOXX future is a large-tick one ( [START_REF] Eisler | The price impact of order book events: market orders, limit orders and cancellations[END_REF]). As a consequence, an upward move in the DAX price (P + D ), while signaling that the market latent price has moved slightly upwards, is not su cient to move the EURO STOXX price by a full tick. However, this move can be perceived in the EURO STOXX through the L b X and C a X Øows that are increasing. Indeed, as already mentioned in Section 4.2, it is well known ([RR10]) that there are more limit orders far from the latent price. The latent price went up, so it is now closer to the best ask, and hence the Øow of the limit (resp. cancel) orders on the best bid (resp. ask) is increasing.

In the opposite direction, a change in EURO STOXX price is perceived as "large" and triggers price changes in the same direction on the DAX. Interestingly, we can also note that changes in the latent price on the EURO STOXX triggers price movements on the DAX. For instance, a shift of liquidity at the bid, namely an increase of the arrival Øow of limit orders at the bid, that signals that the latent price has moved upwards, has a direct e ect on upward price moves on the DAX. This can be seen from the interactions T a

X ! P + D , L b X ! P + D and C a X ! P + D .
We can summarize our results by saying that price changes and liquidity changes on the 6. Conclusion and prospects DAX mainly inØuence liquidity (latent price) on the EURO STOXX, while price changes and liquidity changes on the EURO STOXX tend to trigger price moves on the DAX.

Finally, let us note that the above e ects are even more pronounced when we estimate the interaction matrices with a smaller H . In particular the e ects of DAX price movements on T, L,C on the EURO STOXX become more relevant compared with those on prices. At the same time, while the e ect of EURO STOXX price moves on DAX's ones is still strong, the e ect of liquidity movements on DAX price movements is comparatively stronger with smaller H . This suggests that these e ects are mainly localized at short time scales, while the P + ! P + ones have much slower decay in time.

Bobl -Bund

We perform the same analysis on the asset pair Bobl-Bund futures. Here both assets are large tick assets, however the Bund is much more actively traded than the Bobl in the sense that all the order Øows are of higher intensity. The cross-asset submatrices are depicted in Figure V.11. As in the previous case, we remark that the elements P + L ! P + M and P + M ! P + L reØect the strong correlation observed between the two assets. Price changes in the Bund have also a noticeable e ect on limit/cancel order Øows in the Bobl, while price changes in the Bobl have little to no e ect on the Bund except for the mentioned P + M ! P + L interaction. At the same time, T a , L a ,C a events on the Bobl impact prices on the Bund, while the corresponding event on the Bund have little e ect.

Comparing this with the case of the DAX-EURO STOXX pair, we can liken the e ect of the Bund on the Bobl to that of the DAX over the EURO STOXX and vice-versa. We argue that the di erence in trading frequency between the Bobl and Bund contracts has a similar e ect of that of a di erent tick size that we observed in the previous case. As before, we have an asset, the Bund, which is more "reactive" (the limit/cancel order Øows are higher than those of the Bobl) than the Bobl, thus a price change of the Bund indicating a change of the latent price impacts the limit/cancel Øows of the Bobl. In the previous case, the higher "reactivity" of the DAX was due to its smaller tick size.

Conclusion and prospects

In the context of Hawkes processes, the estimation of the matrix kernel norms is essential, as it gives a clear overview of the dependencies involved in the underlying dynamics. In the context of high-frequency AEnancial time-series non-parametric estimation of the matrix kernel norms has already shown to be very fruitful ( [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF][START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]), since it provides a very rich summary of the system interactions, and it can thus be a valuable tool in understanding a system where many di erent types of events are present. However, its estimation is a computationally demanding process since these estimations are computed from a non-parametric pre-estimation of the kernels themselves, i.e., their entire shape and not only their norm. The resulting complexity prevents the estimations from being performed when the dataset is too heavy or (more important) when the dimension of the Hawkes process (i.e., the number of considered di erent event types) is too large. In this work, we presented the newly developed NPHC algorithm ([ABG + 17]) that allows to directly estimate non-parametrically the kernel norms matrix of a multidimensional Hawkes process, that is without going through the kernel shapes pre-estimation step. As of today, it is the only direct non-parametric estimation procedure available in the academic literature. This method can be seen as a Generalized Method of Moments (GMM) that relies on secondorder and third-order integrated cumulants. This paper shows that this method successfully reveals the various dynamics between the di erent (AErst level) order Øows involved in order books. In a context of a single-asset 8-dimensional Hawkes process, we have shown (as a "sanity check") that it is able to reproduce former results obtained using "indirect" methods. Moreover, the so-obtained gain in complexity allowed us to run a much more detailed analysis (increasing the dimension to 12), separating the di erent types of events that lead to a midprice move. This in turn allowed us to have a very precise picture of the high frequency order book dynamics, revealing, for instance, the di erent interactions that lead to the highfrequency price mean reversion or those between liquidity takers and liquidity makers as well as the inØuence of the tick-size of these dynamics. Not the least, through the analysis of the matrix ™ we also detected the signature of meta-orders. We have also successfully used the NPHC algorithm in a multi-asset 16-dimensional framework. It allowed us to unveil very precisely the high-frequency joint dynamics of two assets that share exposure to the same risk factors but that have di erent characteristics (e.g., di erent tick sizes or di erent degrees of reactivity). It is noteworthy that our methodology can e ciently highlight these types of dynamics, especially since cross-asset e ects are second order e ects compared to mono-asset's.

We conclude by noting that our study left out some relevant information such as the volume of the orders and the size of the jumps in the mid-price. This will be the objective of future works. Moreover, within the methodology presented in this paper, an analysis of baskets of assets (with more than two assets) as well as multi-agent high-frequency interactions are currently under progress.

c W N (µ) = 1 N P N i =1 b g i (µ) b g i (µ) >
, and the objective to minimize is then

√ 1 N N X i =1 b g i (µ) ! °c W N (µ 1 ) ¢ °1 √ 1 N N X i =1 b g i (µ) ! , (20)
where µ 1 is a constant vector. Instead of multiplying by the inverse weighting matrix, we have decided to divide by the sum of its eigenvalues, which is easily computable: 2 for the second. We compute the previous terms with R 1 = 0. All together, the objective function to minimize is APPENDIX A

Tr( c W N (µ)) = 1 N N X i =1 Tr( b g i (µ) b g i (µ) > ) = 1 N N X i =1 Tr( b g i (µ) > b g i (µ)) = 1 N N X i =1 || b g i (µ)|| 2 2 In our case, b g (R) = h vec[ c K c °K c (R)], vec[ b C °C (R)] i > 2 R 2d 2 .
1 k c K c k 2 2 kK c (R) °c K c k 2 2 + 1 k b C k 2 2 kC (R) °b C k 2 2 , (21)

Résumé des contributions

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = E ª [`(µ, ª)],
où f est un terme d'attache dépendant implicitement des données observées,h est un terme de régularisation qui impose une structure à la solution et ª est une variable aléatoire. Typiquement, f est une fonction di érentiable avec un gradient Lipschitz, alors que h pourrait ne pas être lisse -des exemples typiques incluent une pénalité induisant une pénalité -comme la pénalisation `1.

Les algorithmes d'optimisation du premier ordre sont tous des variations de la Descente de Gradient (GD), dont l'origine remonte à Cauchy [START_REF] Cauchy | Méthode générale pour la résolution des systemes d'équations simultanées[END_REF]. A partir d'un point initial µ 0 , cet algorithme minimise une fonction di érentiable f en appliquant itérativement la mise à jour suivante

µ t +1 = µ t °¥t r f (µ t ).
(1) où r f (µ) représente le gradient de f évalué à µ et (¥ t ) est une séquence de tailles de pas. Les algorithmes de descente de gradient stochastique (SGD) se concentrent sur le cas où r f prend beaucoup de temps à calculer, voire est incalculable. En remarquant que r f (µ) s'écrit comme une moyenne, une idée est d'approximer le gradient dans l'étape de mise à jour (1) avec une méthode Monte Carlo par chaîne de Markov [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF]. Par exemple, le remplacement du gradient exact r f (µ) par son estimation MCMC a permis de faire un grand pas en avant dans l'entraînement des modèles graphiques non dirigés [START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF] et des machines Boltzmann restreintes [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. Cette première forme de descente de gradient stochastique est appelée Divergence Contrastive dans les contextes mentionnés.

Introduction

Algorithmes SGD pour une distribution uniforme La plupart des problèmes d'apprentissage statistique de la forme (1) font intervenir comme fonction d'attache aux données f une moyenne sur des points observés, en vertu du principe empirique de minimisation des risques [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. Plus précisément, la fonction objectif écrit

min µ2R d F (µ) = f (µ) + h(µ) with f (µ) = 1 n n X i =1 f i (µ),
où n est le nombre d'observations, et f i est la perte associée à l'observation i th . Dans ce cas, au lieu d'exécuter MCMC pour approcher r f , on échantillonne uniformément un entier aléatoire i entre 1 et r f (µ) et remplace r f (µ) par r f i (µ) dans l'étape de mise à jour (1). Dans la conAEguration à grande échelle, le calcul de r f (µ) à chaque étape de mise à jour représente le goulot d'étranglement de l'algorithme de minimisation, et SGD permet de diminuer le temps de calcul. En supposant que le calcul de chaque r f i (µ) coûte 1, le calcul du gradient complet r f (µ) coûte n, ce qui signiAEe que l'étape de mise à jour de SGD est n fois plus rapide que celle de GD. ) avec AE > 0 est appelé taux de convergence linéaire puisque la diminution de l'erreur après une itération est au pire linéaire. De même, les taux de convergence peuvent être formulés comme la complexité totale pour atteindre une précision AExe, c'est-à-dire le nombre d'itérations après lequel la di érence

E f (µ t ) °f (µ t ) °f (µ §

Processus ponctuels

Le processus ponctuel est un outil mathématique utile pour décrire les phénomènes qui se produisent à des endroits et/ou à des moments aléatoires. Un processus de points est un élément aléatoire dont la valeur est une liste de points sur un ensemble S. Nous présentons ici les résultats utiles lorsque l'ensemble S est l'intervalle [0, T ), et les points sont des événements datés; ce cas spécial est parfois appelé processus de point temporel. Le livre [START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] est considéré comme la référence principale sur la théorie des processus ponctuels. Deux exemples de processus ponctuels temporels sont traités dans cette thèse. Le premier est le processus ponctuel derrière le modèle proportionnel des risques de Cox : sa fonction d'intensité conditionnelle permet de déAEnir le hazard ratio, une quantité fondamentale dans la littérature d'analyse de survie, voir [START_REF] Andersen | Statistical models based on counting processes[END_REF]. Le modèle de régression de Cox relie la durée avant un événement appelé échec à certaines covariables. Ce modèle peut être reformulé dans le cadre de processus ponctuels [START_REF] Andersen | Statistical models based on counting processes[END_REF]. Le second est le processus de Hawkes qui modélise comment les événements passés augmentent la probabilité d'événements futurs. Sa version multivariée permet d'encoder une notion de causalité entre les di érents noeuds. Nous présentons ci-dessous le modèle des risques proportionnels de Cox et les processus de Hawkes dans la partie II.

Modèle des risques proportionnels de Cox

L'analyse de survie étudie la durée qui précède l'arrivée d'un évènement particulier, tel que la mort dans les organismes biologiques et les défaillances dans les systèmes mécaniques, et est maintenant répandue dans une variété de domaines comme la biométrie, l'économétrie et l'assurance. La variable que nous étudions est le temps d'attente jusqu'à ce qu'un événement Introduction bien déAEni se produise, et l'objectif principal de l'analyse de survie est de lier les covariables, ou caractéristiques, d'un patient à son temps de survie T .

Suivant la théorie des processus ponctuels, nous déAEnissons l'intensité comme la probabilité conditionnée qu'un patient meurt immédiatement après t , étant donné qu'il était vivant avant t :

∏(t ) = lim h!0 P(t ∑ T ∑ T ∑ t + h|t ∑ T ) h .
L'approche la plus populaire, pour certaines raisons expliquées ci-dessous, est le modèle des risques proportionnels de Cox [START_REF] David | Regression models and life tables (with discussion)[END_REF]. Le modèle de Cox prend une forme semiparamétrique pour le ratio de risque au temps t pour le patient i , dont les caractéristiques sont codées dans le vecteur x i 2 R d :

∏ i (t ) = ∏ 0 (t ) exp(x > i µ),
où ∏ 0 (t ) est un ratio de risque de base, qui peut être considéré comme le ratio de risque d'un patient dont les covariables sont x = 0. Une approche d'estimation considère ∏ 0 comme une nuisance et estime seulement µ en maximisant une vraisemblance partielle [START_REF] David | Regression models and life tables (with discussion)[END_REF]. Cette façon d'estimer convient aux études cliniques où les médecins ne s'intéressent qu'aux e ets relatifs des covariables codées en x sur le ratio de risque. Pour ce faire, on peut calculer le rapport des rapports de risque de deux patients di érents :

∏ i (t ) ∏ j (t ) = exp((x i °x j ) > µ)
Pour cette raison, on dit que le modèle de Cox est un modèle de risques proportionnels.

SVRG au delà de la Minimisation du Risque Empirique

Les données utilisées en analyse de survie de la forme (y Notre algorithme de minimisation est doublement stochastique dans le sens où les étapes de gradient sont faites en utilisant la descente stochastique de gradient (SGD) avec réduction de variance, et les espérances internes sont approximées par un algorithme de la chaîne de Monte Carlo Markov (MCMC). Nous dérivons des conditions sur le nombre d'itérations MCMC garantissant la convergence, et obtenons un taux de convergence linéaire sous forte convexité et un taux sublinéaire sans cette hypothèse.

i , x i , x i , ± i ) i =1 i =1 n pat contiennent

Part II: Découvrir la causalité de Hawkes sans paramétrage

Dans les Chapitres III et IV, nous étudions deux méthodes permettant de retrouver les relations de causalité à partir d'un processus ponctuel multivarié. Nous développons une approche par chapitre.

Processus de Hawkes

AAEn de modéliser la dynamique commune de plusieurs processus ponctuels (par exemple l'horodatage des messages envoyés par di érents utilisateurs d'un réseau social), nous allons considérer le modèle multi-dimensionnel de Hawkes, introduit en 1971 dans [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF] et [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF], avec des inØuences croisées entre les di érents processus. Par déAEnition, une famille de d processus ponctuels est un processus de Hawkes multi-dimensionnel si les intensités de toutes ses composantes s'écrivent comme des régressions linéaires sur le passé des processus d : qui compte le nombre d'événements de i dont l'ancêtre direct est un événement de j , on sait d'après [START_REF] Bacry | Hawkes processes in AEnance[END_REF] que : Dans la littérature, il existe deux principales classes de procédures d'estimation pour les noyaux de Hawkes : la paramétrique et la non paramétrique. La première suppose une paramétrisation des noyaux de Hawkes, la plus courante suppose que les noyaux sont en décomposition exponentielle, et estime le paramètre via la maximisation de la log-vraisemblance de Hawkes, voir par exemple [START_REF] Bacry | A generalization error bound for sparse and low-rank multivariate hawkes processes[END_REF] ou [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF].La seconde est basée soit sur la résolution numérique des équations de Wiener-Hopf qui relie les noyaux de Hawkes à sa structure de corrélation [?] (ou de manière équivalente sur l'approximation du processus de Hawkes en tant que modèle autorégressif et la résolution des équations de Yule-Walker [START_REF] Eichler | Graphical modeling for multivariate hawkes processes with nonparametric link functions[END_REF]), soit sur une méthode de moments via la minimisation de la fonction de contraste déAEnie dans [START_REF] Reynaud-Bouret | Goodness-of-AEt tests and nonparametric adaptive estimation for spike train analysis[END_REF].

∏ i t = µ i + D X k=1 Z t 0 ¡ i j (t
E[d N i √ j t ] = g i j E[d N j t ] = g i j § j d t,
Dans les Chapitres III et IV, nous proposons deux méthodes d'estimation non-paramétrique permettant d'estimer les intégrales des noyaux de Hawkes à l'aide des intégrales des moments du processus.

Pour toutes les procédures d'estimation mentionnées ci-dessus, y compris la nôtre, nous avons besoin de la condition de stabilité suivante aAEn que le processus admette une version avec une intensité stationnaire : Assumption 1. La norme spectrale of G = [g i j ] satisfait ||G|| < 1.

Approche par Méthode des Moments Généralisée

Un travail récent [JHR15] a prouvé que les cumulants intégrés des processus de Hawkes peuvent être exprimés en fonctions de G = [g i j ], et a fourni la méthode constructive pour obtenir ces expressions. La première approche que nous avons développée dans cette partie est une méthode d'appariement des moments sur les cumulants intégrés de deuxième et de troisième ordre du processus. À cette AEn, nous avons conçu des estimateurs cohérents des premier, deuxième et troisième cumulants intégrés du processus de Hawkes. Leurs contrepar-2. Part II: Découvrir la causalité de Hawkes sans paramétrage ties théoriques sont des polynômes de R = (I °G) °1, comme indiqué dans [START_REF] Jovanović | Cumulants of hawkes point processes[END_REF] :

§ i = d X m=1 R i m µ m C i j = d X m=1 § m R i m R j m K i jk = d X m=1 (R i m R j m C km + R i m C j m R km +C i m R j m R km °2 § m R i m R j m R km ).
Une

Modèle 12-dimensionnel du carnet d'ordres d'un actif

Comme première application de la procédure décrite au chapitre III, nous considérons le processus ponctuel à 12 dimensions suivant, une extension naturelle du processus ponctuel à 8 dimensions introduit dans [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] : • T + (T °): mouvement du mid-price à la hausse (baisse) due à un ordre de marché.

N t = (T + t ,
• L + (L °): mouvement du mid-price à la hausse (baisse) due à un ordre limite.

• C + (C °): mouvement du mid-price à la hausse (baisse) due à un ordre d'annulation.

• T a (T b ): ordre de marché à l'ask (au bid) sans modiAEcation du mid-price.

• L a (L b ): ordre limite à l'ask (au bid) sans modiAEcation du mid-price.

• C a (C b ): ordre d'annulation à l'ask (au bid) sans modiAEcation du mid-price.

Nous utilisons ensuite l'interprétation causale des processus de Hawkes pour interpréter notre solution comme une mesure de la causalité entre les événements. Cette application de la méthode à ce nouveau modèle a révélé les di érentes interactions qui conduisent au retour à la moyenne des prix à haute fréquence, et celles entre les preneurs de liquidité et les faiseurs de liquidité.

Par exemple, on observe les e ets des événements T + sur d'autres événements sur la

Modèle 16-dimensionnel du carnet d'ordres de deux actifs

La méthode d'estimation non paramétrique introduite au chapitre III permet une estimation rapide pour une méthodologie non paramétrique. Nous passons ensuite à l'échelle le modèle aAEn de tenir compte des événements sur deux actifs simultanément et de dévoiler une structure précise de la dynamique des actifs croisés à haute fréquence. Nous considérons un modèle à 16 dimensions, composé de deux modèles à 8 dimensions de la forme. Abstract : This thesis is divided into three parts. The first focuses on a new optimization algorithm we have developed. It allows to estimate the parameter vector of the Cox regression when the number of observations is very important. Our algorithm is based on the SVRG algorithm and uses a MCMC method to approximate the descent direction. We have proved convergence rates for our algorithm and have shown its numerical performance on simulated and real world data sets. The second part shows that the Hawkes causality can be estimated in a non-parametric way from the integrated cumulants of the multivariate point process. We have developed two methods for estimating the integrals of the kernels of the Hawkes process, without making any hypothesis about the shape of these ker-nels. Our methods are faster and more robust, with respect to the shape of the kernel, compared to the state-of-the-art. We have demonstrated the statistical consistency of the first method, and have shown that the second method can be reduced to a convex optimization problem. The last part highlights the dynamics of the order book thanks to the first nonparametric estimation method introduced in the previous section. We used EUREX futures data, defined new order book models and applied the estimation method on these point processes. The results obtained are very satisfactory and consistent with an econometric analysis, and prove that the method that we have developed makes it possible to extract a structure from data as complex as those resulting from high-frequency finance.

N

1.

  Part I: Large-scale Cox model typical examples include sparsity inducing penalty -such as the `1 penalization.

  Introduction which the di erence E f (µ t ) °f (µ § ) becomes smaller than ≤ > 0 multiplied by the complexity per iteration. The algorithm Gradient Descent will reach the accuracy ≤ after O °∑ log 1 ≤ ¢ iterations resulting in a O °nd∑ log 1 ≤ ¢ complexity, while Stochastic Gradient Descent reachs such accuracy after O °∑ ≤

  to reach accuracy ≤) than Gradient Descent since the complexity per iteration of those algorithms is O(d ) versus O(nd) for Gradient Descent. Those typically typically enjoy complexity of the form O °(n + ∑)d log 1 ≤ ¢ in the strongly-convex case, see [SLRB17, JZ13, DBLJ14, SSZ13].

Result 1 .

 1 Under Assumptions 1 and 2, the sequence of estimators deAEned by the minimization of L T (R) converges in probability to the true value G:

  is a norm that provides a particular structure to the solution. Every matrix G satisfying C = (I °G) °1L(I °G> ) °1 equals I °L1/2 MC °1/2 with M an orthogonal matrix. Instead of the previous problem, we now focus on its convex relaxation, we split the variables G and M, and solve the problem with the Alternating Direction Method of Multipliers algorithm, see [GM75] and [GM76]: min G,M f (G) + B (M) + B (G) + R d £d + (G) s.t. G = I °L1/2 M C °1/2 , where B (resp. B) is the open (resp. closed) unit ball w.r.t. the spectral norm. The closed unit ball w.r.t. the spectral norm is indeed the convex hull of the orthogonal group.

  Figure .1: Kernel norm matrix G estimated for the DAX future with H = 1s.

  which basically is the number of events in each Borel subset B 2 B. The mean measure M of a point process ª is a measure on S that assigns to every B 2 B the expected number of events of ª in B , i.e., M (B ) := E[N (B )] for all B 2 B. For inhomogeneous Poisson process, M (B ) = R B ∏(x)d x, where the intensity function ∏(x) yields a positive measurable function on S. Intuitively speaking, ∏(x)d x is the expected number of events in the inAEnitesimal d x. For the most common type of point process, a homogeneous Poisson process, ∏(x) = ∏ and M (B ) = ∏|B |, where | • | is the Lebesgue measure on (S, B). More generally, we deAEne Cox point processes -also known as doubly stochastic Poisson processes -as a generalization of Poisson processes where the intensity ∏(x) is itself a stationary stochastic process. Then, conditional on ∏, the doubly stochastic Poisson process is simply an inhomogenous Poisson process with intensity ∏(x).

) 4 .

 4 Theoretical guarantees for the iterations t belonging to the k-th phase, where N k is the number of iterations of the Markov chain used for the computation of b r f i t (µ t ) during phase k (see line 5 of Algorithm 6), and where C 1 and C 2 are positive constants. Let us point out that Proposition 1 below gives a su cient condition for Assumption 2 to hold.

Figure

  Figure II.1: Convergence of Cocktail and L-BFGS-B on Lymphoma dataset. Top: the starting point is µ 0 = 0 2 R d . Bottom: the starting point is µ 0 = b µ (l ) (solution to the same objective with a slightly larger ∏). This illustrates the fact that Cocktail cannot minimize directly a single objective (with a AExed ∏) and requires to compute the full path of solution to converge.

  Figure II.1, where the convergence of Cocktail and L-BFGS-B algorithms are compared for two starting points µ 0 . Even when the starting point is set to the previous minimizer (second case in Figure II.1, cocktail's convergence is slower than the one of L-BFGS-B. As a consequence, we decided that no fair comparison could be conducted with cocktail and coxnet algorithms.

Figure

  Figure II.2: Distance to optimum of all algorithms on NKI70 (left) and Lymphoma (right) with ridge penalization (AE = 0 and ∏ = 1/ p n)

pn

  and illustrate our results in Figure II.5. Low lasso. We take AE = 1 and ∏ = 1/n and illustrate our results in Figure II.4. High ridge. We take AE = 0 and ∏ = 1/ p n and illustrate our results in Figures II.2 and Figures II.3. Low ridge. We take AE = 0 and ∏ = 1/n and illustrate our results in Figure II.6.

Figure

  Figure II.4: Distance to optimum of all algorithms on NKI70, Lymphoma, Luminal and on the simulated dataset (respectively from top to bottom) for Low-ridge penalization

-

  increasing: no memory process and patient's health is worsening This method enables us to simulate n failures times T 1 , T 2 ,...,T n . Then, we simulate C 1 ,C 2 ,...,Cn with exponential distribution. This AEnally gives us a set of observed times (yi ) n i =1 = (T i ^Ci ) n i =1 and a set of censoring indicators (± i ) n i =1 = ( {T i ^Ci } ) n i =1 .10 Mini-batch sizingThe mini-batch sizing question is essential since it is a natural trade-o between computing time and precision. We know that computing r f i (µ) needs the computation of |R i | 2 {1, . . . , n pat } inner products. One proves easily that computing a mini-batch (1/n mb )r f B (µ)where B is the set of n mb index randomly picked -only needs max i 2B |R i | inner products. A simple probability exercise gives us a key insight about the mini-batch size. Let's assume that censoring is uniform over the set {1, 2, . . . , n pat } meaning that |R i | = ci with c > 1. Then, we denote u 1 , u 2 ,...,u n mb ª U [n] the indices independently sampled to compute the mini-batch i.e. B = {u i } n mb i =1 . Now we study the c.d.f. of max 1∑i ∑n mb u i : for k 2 {1, 2, . . . , n},

i

  2B |R i | and |B| from being too large. This is why we used n mb = 0.1n or n mb = 0.01n, depending of the size n of the dataset. 54 10. Mini-batch sizing Part II

  b § i , b C i j , b K i i j fromEqs. (11, 12, 13) 4: Design f L (R) using the computed estimators. 5: Minimize numerically f L (R) so as to obtain b R 6: Return b G = I d °b R °1.

  ) writes C 1/2 OL °1/2 , with L = diag( §) and O an orthogonal matrix. Our starting point is then simply chosen by setting O = I d in the previous formula, leading to nice convergence results. Even though our main concern is to retrieve the matrix G, let us notice we can also obtain an estimation of the baseline intensities' from Eq. (3), which leads to b µ = b R °1 b §. An e cient implementation of this algorithm with TensorFlow, see [AAB + 16], is available on GitHub: https://github.com/achab/nphc.

  can be seen as the weighted squared Frobenius norm of b g T (R).Moreover, when T ! +1, one has b under the conditions of the following theorem, where P ! stands for convergence in probability.

Figure III. 1 :

 1 Figure III.1: The three di erent kernels used to simulate the datasets.

Figure III. 2 :

 2 Figure III.2: On Exp100 dataset, estimated bG with ADM4 (left), with NPHC (middle) and the ground-truth matrix G (right). Both ADM4 and NPHC estimates recover the three blocks. However, ADM4 overestimates the integrals on two of the three blocks, while NPHC gives the same value on each blocks.

Figure

  Figure III.3: Estimated b G via NPHC on DAX order book data.

  Considering a block-wise weighting matrix, one block for c K c °K c (R) and the other for b C °C (R), the sum of the eigenvalues of the AErst block becomes k c K c °K c (R)k 2 2 , and k b C °C (R)k 2

Theorem 3 .

 3 Suppose that (N t ) is observed on R + , cW T P °! W , and 1. W is positive semi-deAEnite and W g 0 (µ) = 0 if and only if

Remark 1 .

 1 and || f || 1 = O(1). In practice, we use a constant sequence of weighting matrices: c W T = I d . Proof. Proceed by verifying the hypotheses of Theorem 2.1 from [NM94]. Condition 2.1(i ) follows by (i ) and by Q 0

  µ)| P °! 0, we should now prove that sup µ2£ k b g T (µ) °g0 (µ)k P °! 0. By £ compact, it is su cient to prove that k b

  (skewness density) and M i jkl u,v,w (fourth cumulant density) are positive as polynoms of integrals of √ ab • with positive coe cients. The integrals of the singular parts are positive as well.

Figure IV. 1 :

 1 Figure IV.1: From the left to the right: solution of Problem I, solution of Problem II, solution of Problem III, and the ground-truth matrix G. Only the solution to Problem I recovers the three blocks. The two other problems outputed symmetric matrices, while the ground-truth matrix is not.

Figure IV. 2 :

 2 Figure IV.2: Solutions of Problem I for the multivariate point process A t (left) and B t (right). We observe a strong self-excitation. These solutions are consistent with the estimated kernel norm matrix in [RBL17].

Figure IV. 3 :

 3 Figure IV.3: Solution of Problem I for the multivariate point process Ct . This solution is consistent with estimtates in lower dimension.

  Figure V.1: The two di erent kernels used to simulate the datasets.

Figure V. 3 :

 3 Figure V.3: Kernel norm matrix G estimated with the NPHC method for the DAX future (left) and with the Wiener-Hopf method of [BM16] (right) when the 8-dimensional model described in Section 4.2 is considered.

Figure V. 4 :

 4 Figure V.4: Kernel norm matrix G for the DAX future estimated (using NPHC) at three di erent times: between 08:00 and 10:00 (left), between 12:00 and 14:00 (middle) and between 16:00 and 18:00 (right).

Figure V. 5 :

 5 Figure V.5: Estimation of the baseline intensities of each event type within a trading day for the DAX future using 15 min slots. Left panel: Empirical intraday pattern measured using market, limit and cancel orders that do not move the price. Right panel: µ values estimated using the NPHC method. All quantities are expressed in s °1.

110 4 .

 4 Figure V.6: Kernel norm matrix G estimated with the NPHC method for the DAX future with H = 1s (left) and H = 500s (right).

Figure V. 7 :

 7 Figure V.7: Kernel norm matrix G estimated with the NPHC method for the Bund future with H = 1s (left) and H = 500s (right).

Figure

  Figure V.8: ™ matrix of eq. (6) estimated with the NPHC method for the DAX future (left) and the Bund future (right) with H = 500s.

Figure V. 9 :

 9 Figure V.9: Hawkes kernel norm matrix obtained when the DAX and EURO STOXX futures are considered simultaneously in a 16D model. DAX events are denoted with the D subscript, EURO STOXX ones with the X subscript.

Figure V. 10 :

 10 Figure V.10: Submatrices of the Kernel norm matrix G corresponding to the e ect of DAX events on EUROSTOXX STOXX events (left) and vice versa (right). These two submatrices correspond to the ones lying on the antidiagonal on the Figure V.9

  Assuming the associated weighting matrix is block-wise, one block for c K c °K c (R) and the other for b C °C (R), the sum of the eigenvalues of the AErst block becomes k c K c °K c (R)k 2 2 , and k b C °C (R)k 2

  Figure V.11: Submatrices of the Kernel norm matrix G corresponding to the e ect of Bund (L) events on Bobl (M) events (left) and vice-versa (right).

)

  devient plus petit que ≤ > 0 multiplié par la complexité par itération. L'algorithme Descente de Gradient atteindra la précision ≤ après O °∑ log 1 ≤ ¢ itérations résultant en une complexité de O °nd∑ log 1 ≤ ¢ , tandis que la Descente de Gradient Stochastique atteint une telle précision après O °∑ ≤ ¢ itérations et donc une complexité en O ≥ d ∑ ≤ ¥ . Récemment, di érents travaux ont amélioré la descente du gradient stochastique en utilisant les techniques de réduction de variance des méthodes de Monte Carlo. L'idée est d'ajouter un terme de contrôle à la direction de descente pour améliorer le compromis biais-variance dans l'approximation du gradient réel r f (µ). Ces variantes bénéAEcient également de taux de convergence linéaire, puis de complexités plus petites (pour atteindre la précision ≤) que la descente de gradient puisque la complexité par itération de ces algorithmes est de O(d ) contre O(nd) pour Gradient Descent. Celles-ci atteignent typiquement une complexité de la forme O °(n + ∑)d log 1 ≤ ¢ dans le cas fortement convexe, voir [SLRB17, JZ13, DBLJ14, ?]. 1. Part I: Modèle de Cox à grande échelle

∏t

  = ∏ > 0. Notez que les processus ponctuels temporels peuvent aussi être caractérisés par la distribution des temps d'intervalle i.e. la durée entre deux événements consécutifs. Nous rappelons que la distribution des temps d'intervalle d'un processus de Poisson avec intensité ∏ est une distribution exponentielle du paramètre ∏. Voir la page 41 de[START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] pour quatre façons équivalentes de déAEnir un processus ponctuel temporel.

  Cependant, maximiser cette vraisemblance partielle est un problème di cile lorsqu'il s'agit de données à grande échelle (c'est-à-dire un grand nombre d'observations n) et à haute dimension (c'est-à-dire un grand d ). Pour s'attaquer à la dimensionnalité élevée, des approches pénalisées et parcimonieuses ont été envisagées dans la littérature[START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] [T + 97][START_REF] Goeman | L1 penalized estimation in the cox proportional hazards model[END_REF]. Le problème est maintenant de minimiser l'opposé du logarithme de la vraisemblance partielle f (µ) = °`(µµ) avec une pénalisation h(µ) qui fait que le prédicteur µ devient parcimonieux et sélectionne des variables. Nous discuterons de cette approche et des di érents modèles dans le chapitre II. Il n'existe cependant pas encore d'approches pour répondre au problème à grande échelle.

  , pour chaque individu i = 1, . . . , n pat , un vecteur de caractéristiques x i 2 R d , un temps observé y i 2 R + qui correspond au temps de l'échec si ± i = 1 ou à un temps censuré à droite si ± i = 0. Si D = {i : ± i = 1} est l'ensemble des patients pour lesquels un temps d'échec est observé, si n = |D| est le nombre total de temps d'échec, et si R i = { j... : y j ∏ y i } est l'indice des individus 8 2. Part II: Découvrir la causalité de Hawkes sans paramétrage toujours à risque au moment où y i , l'opposé du logarithme de la vraisemblance partielle de Cox s'écrit : µ 2 R d . Chaque gradient de la probabilité logarithmique négative s'écrit alors comme deux espérances imbriquées : l'une d'une distribution uniforme sur D, l'autre sur une distribution de Gibbs, voir le chapitre II pour plus de détails.

  § i correspond à l'intensité moyenne, satisfaisant E[d i t ] = § i d t.Cependant, dans la pratique, les noyaux de Hawkes ne sont pas directement mesurables à partir des données et ces mesures de causalité entre les di érents types d'événements sont donc inaccessibles.

  ∏ 0 où f (G) est une norme qui donne une certaine structure à la solution. Toute matrice G satisfaisant C = (I °G) °1L(I °G> ) °1 s'écrit I °L1/2 MC °1/2 avec M une matrice orthogonale. Au lieu d'étudier le problème précédent, nous nous focalisons sur sa relaxation convexe, we séparons les variables G et M, et résolvons le problème avec l'algorithme Alternating Direction Method of Multipliers, voir [GM75] et [GM76]: min G,M f (G) + B (M) + B (G) + R d £d + (G) s.t. G = I °L1/2 M C °1/2 , où B (resp. B) est la boule unitaire ouverte (resp. fermée) par rapport à la norme spectrale. La boule unitaire fermée par rapport à la norme spectrale est en e et l'enveloppe convexe du groupe orthogonal. Contrairement au problème d'optimisation du chapitre précédent, le problème qui vient d'être énoncé est convexe. Nous testons cette procédure sur des simulations numériques de divers noyaux Hawkes et des données du carnet d'ordres réels, et nous montrons comment le critère f a ecte les matrices que nous récupérons. 3 Partie III: Capter la dynamique d'un carnet d'ordres à l'aide de processus de Hawkes Le chapitre V s'intéresse à l'estimation des intégrales des noyaux de Hawkes sur des données AEnancières, à l'aide de la méthode d'estimation introduite dans le chapitre III. Cela nous a permis d'avoir une image très précise de la dynamique du carnet d'ordres à haute fréquence. Nous avons utilisé les événements du carnet de commandes associés à 4 actifs très liquides de la bourse EUREX, à savoir DAX, EURO STOXX, Bund et les contrats à terme Bobl.

Figure A. 1 (

 1 Figure A.1: Matrice d'intégrales des noyaux G estimée pour le DAX avec H = 1s.

Figure A. 2 :

 2 Figure A.2: Sous-matrice d'intégrales des noyaux G correspondant à l'e et des évènements du DAX sur ceux de l'EUROSTOXX STOXX (gauche) et vice versa (droite).

  • NKI70 contains survival data for 144 breast cancer patients, 5 clinical covariates and the expressions from 70 gene signatures, see [VDVHVV + 02].

• Luminal contains survival data for 277 patients with breast cancer who received the adjuvant tamoxifen, with 44,928 expressions measurements, see [LHKD + 07]. • Lymphoma contains 7399 gene expressions data for 240 lymphoma patients. The data was originally published in [AED + 00].

  causality if forecasting future values of Y is more successful while taking X past values into account. In[START_REF] Eichler | Graphical modeling for multivariate hawkes processes with nonparametric link functions[END_REF], it is shown that for N

t a multivariate Hawkes process, N j t does not Granger-cause N i t w.r.t N t if and only if ¡ i j

Table III .

 III 1: Complexity of state-of-the-art methods. NPHC's complexity is very low , especially in the regime n ¿ d .

	Method	Total complexity
	ODE [ZZS13]

Table III

 III

	Method	ODE GC	ADM4 WH	NPHC
	RelErr	0.007 0.15 0.10	0.005 0.001
	MRankCorr 0.33	0.02 0.21	0.34	0.34
	Time (s)	846	768	709	933	20
	Table III.3: Metrics on PLaw10: comparable rank correlation, strong improvement for relative
	error and computing time.
	Method	ODE GC	ADM4 WH	NPHC
	RelErr	0.011 0.09 0.053	0.009 0.0048
	MRankCorr 0.31	0.26 0.24	0.34	0.33
	Time (s)	870	781	717	946	18

.2: Metrics on Rect10: comparable rank correlation, strong improvement for relative error and computing time.

Table III

 III

	.4: Metrics on Exp100: comparable rank correlation, strong improvement for relative
	error and computing time.
	Method	ODE	GC	ADM4 NPHC
	RelErr	0.092 0.112	0.079	0.008
	MRankCorr 0.032 0.009 0.049 0.041
	Time (s)	3215	2950	2411	47
	Table III.5: Metrics on MemeTracker: strong improvement in relative error, rank correlation
	and computing time.
	Method	ODE GC	ADM4 NPHC
	RelErr	0.162 0.19	0.092	0.071
	MRankCorr 0.07	0.053 0.081	0.095
	Time (s)	2944 2780 2217	38

  ||•|| 1 , the squared `2-norm ||•|| 2 2 and the nuclear norm ||•|| § . In the rest of the section, we refer as Problem I the minimization problem written in Equation (2) with

Table IV .

 IV 1: The solution of Problem I has the smallest relative error.

	Problem	I	I I	I I I
	RelErr	0.093 0.130 0.131

1 Part I: Modèle de Cox à grande échelle

  Cette thèse cherche à montrer comment certaines méthodes d'optimisation récentes permettent de résoudre des problèmes d'estimation di ciles liés aux modèles d'évènements. Alors que le cadre classique de l'apprentissage supervisé[START_REF] Hastie | Overview of supervised learning[END_REF] traite les observations comme une collection de couples indépendants de covariables et de labels, les modèles d'évènements s'intèressent aux temps d'arrivée d'un évènement et cherchent alors à extraire de l'information de la source de donnée. Ces évènements datés sont ordonnés par la chronologie et ne peuvent dès lors être considérés comme indépendants. Ce simple constat motive l'utilisation d'un outil mathématique particulier appelé processus ponctuel[START_REF] Daley | An introduction to the theory of point processes: volume II: general theory and structure[END_REF] pour apprendre une structure à partir de ces évènements. Nous allons dans un premier temps présenter et motiver les problématiques que nous voulons aborder dans cette thèse. Le fait que tous les patients ne meurent pas lors de l'étude est toujours intéressant d'un point de vue statistique, mais ces données ne peuvent pas être utilisée dans un problème de régression classique pour lequel il faudrait observer l'évènement défaillance pour tous les individus. La di culté a été contournée par D.R. Cox[START_REF] David | Regression models and life tables (with discussion)[END_REF], dans l'un des articles scientiAEques les plus cités de tous les temps[START_REF] Van Noorden | The top 100 papers[END_REF] avec le modèle à risque proportionnel qui permet d'extraire de l'information de données censurées, i.e. de patients pour lesquels le temps de défaillance n'est pas observé. La procédure d'estimation du vecteur de paramètre de la régression, sans aucune hypothèse sur le risque de base considéré comme un paramètre de bruit, a été introduite dans[START_REF] Cox | Partial likelihood[END_REF] Introduction et revient à maximiser la vraisemblance partielle du modèle. Cette procédure permet de gérer e cacement les covariables de grande dimension, ce qui est courant avec les données de biostatistique, en ajoutant un terme de pénalisation au critère à minimiser[START_REF] Goeman | L1 penalized estimation in the cox proportional hazards model[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Cependant, les algorithmes de maximisation de la vraisemblance partielle ne passent pas à l'échelle lorsque le nombre de patients devient très grand, contrairement à la plupart des algorithmes qui ont permis la data révolution. On peut dès lors se poser la question suivante :Quelques années avant le vingtième siècle, le sociologue français Durkheim a rmait déjà que les sociétés humaines sont faites de composantes interconnectées, comme les systèmes biologiques[START_REF] Durkheim | Le suicide: étude de sociologie[END_REF]. Maintenant que notre technologie nous permet même d'être connecté à distance, la notion de réseau concerne de très nombreux domaines : réseaux sociaux, systèmes d'information, marketing, épidémiologie, sécurité nationale et tant d'autres. Une meilleure compréhension de ces larges réseaux et des processus qui s'y passent aurait des applications majeures dans les domaines déjà cités[START_REF] Rodriguez | Structure and Dynamics of Di usion Networks[END_REF]. L'observation des réseaux se réduit souvent à l'enregistrement des instants où les noeuds du réseau envoient un message, achètent un produit ou sont infectés par un virus. Nous observons souvent où et quand mais pas comment et pourquoi les messages sont envoyés via un réseau social. L'obtention de ces données pour plusieurs noeuds du réseau permet de retrouver la dynamique jointe et révéler la structure sous-jacente au système. Une des approches permettant d'estimer l'inØuence entre ces di érentes sources est d'utiliser un processus ponctuel appelé processus de Hawkes[START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF][START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF], dont le taux d'arrivée des évènements dépend des évènements passés. Les processus de Hawkes ont été appliqués avec succès pour modéliser l'inØuence réciproque entre les tremblements de terre de di érentes magnitudes et rapprochés dans le temps[START_REF] Ogata | Statistical models for earthquake occurrences and residual analysis for point processes[END_REF]. Plus précisement, ce processus quantiAEe l'augmentation de la probabilité d'observer de nouveaux tremblements de terre, appelés répliques, après en avoir observé un premier, via l'utilisation de fonctions appelés noyaux. Les processus de Hawkes permettent aussi de mesurer la causalité au sens de Hawkes, qui correspond au number moyen d'évènements de type i engendrés par un évènement de type j . Outre l'exemple originel des tremblements de terre, les deux autres domaines majeurs où sont utilisés les processus de Hawkes sont l'étude des réseaux sociaux [BBH12, ZZS13, ISG13] et l'étude des transactions AEnancières[START_REF] Bacry | Hawkes processes in AEnance[END_REF]. L'estimation habituelle de la causalité au sens de Hawkes demande cependant de faire quelques hypothèses sur la formes des noyaux pour simpliAEer l'algorithme d'inférence[START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF]. Une hypothèse courante est la décroissance monotone des noyaux (exponentiels ou en loi du puissance), ce qui signiAEe que l'impact d'un événement est toujours instantanément maximal, ce qui n'est pas réaliste puisqu'en pratique il peut exister un délai avant l'impact maximal. Ces remarques nous conduisent à soulever la question suivante : 'échelle lorsque le nombre de noeuds est élevé. Dans cette partie, nous ne nous concentrons que sur la première, pour laquelle nous avons prouvé un résultat de consistence. Depuis le travail pionnier de Bowsher[START_REF] Bowsher | Modelling security market events in continuous time: Intensity based, multivariate point process models[END_REF], qui a reconnu la Øexibilité et la simplicité d'utilisation des processus Hawkes aAEn de modéliser la dynamique conjointe des transactions et des changements de prix moyens du NYSE, les processus Hawkes n'ont pas cessé de gagner en popularité dans le domaine de la AEnance à haute fréquence, voir[START_REF] Bacry | Hawkes processes in AEnance[END_REF] pour une revue. En e et, pour prendre en compte des données transactionnelles irréguliers espacées dans le temps, il est naturel de les considérer comme un processus ponctuel. Aussi, dans le domaine AEnancier, de nombreuses caractéristiques résumant les résultats empiriques sont déjà connues. Par exemple, le Øux des transactions est connu pour être auto-corrélé et inter-corrélé avec les mouvements de prix. Ces caractéristiques, appelées faits stylisés, de l'économiste Nicholas Kaldor[START_REF] Kaldor | A model of economic growth[END_REF], faisaient référence à des tendances statistiques qui doivent être prises en compte malgré un possible manque de compréhension microscopique. Ces faits stylisés peuvent être facilement capturés à l'aide de la notion de causalité de Hawkes. La compréhension de la dynamique du carnet d'ordre est l'une des questions centrales des statistiques AEnancières, et les représentations non-paramétriques antérieures des carnets d'ordre à l'aide de processus de Hawkes multivariés étaient de faible dimension en raison de la complexité de leur méthode d'estimation. L'estimation non-paramétrique de la causalité de Hawkes introduite dans la deuxième partie de cette thèse est rapide et robuste à la forme des fonctions noyaux, et il est donc naturel de se demander quel type de fait stylisé peut être découvert à partir des données horodatées du carnet d'ordre.Dans la partie I, nous répondons à la question 4 en introduisant un nouvel algorithme de descente de gradient stochastique appliquée à la minimisation de la vraisemblance partielle de Cox. En e et, la log-vraisemblance partielle de Cox s'écrit comme une somme de sousfonctions dépendant chacune d'une séquence d'observations, séquences de longueur variable, contrairement au cas classique de la minimisation du risque empirique où les sous-fonctions dépendent d'un nombre AExe d'observations, une en général. Les algorithmes classiques de descente de gradient stochastique sont moins e caces dans notre cas. Nous avons adapté l'algorithme SVRG[START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] [XZ14] en rajoutant une nouvelle étape d'échantillonage : chaque sous-fonction est approximée par une méthode de Monte Carlo par chaînes de Markov (MCMC), son calcul exact étant coûteux. Notre algorithme jouit d'un taux de convergence linéaire, une fois que le nombre d'itération de la chaîne de Markov est plus grand qu'une borne inférieure explicite. Nous illustrons la surperformance de notre algorithme sur des jeux de données issus de l'analyse de survie.IntroductionLes réponses à la question 5 se trouvent dans la partie II où nous étudions deux algorithmes d'estimation non-paramétrique de l'estimation de la causalité de Hawkes. Ces deux méthodes se basent sur le calcul des cumulants intégrés du processus de Hawkes multivarié et tirent parti des relations polynomiales entre ces cumulants intégrés et la matrice de causalité de Hawkes. La première approche repose sur la correspondance entre l'écriture théorique et le calcul empirique des cumulants du deuxième et troisième ordre. Cela se fait par la minimisation de la norme quadratique de la di érence entre les deux termes, ce qui peut être vu comme un cas de Méthode des Moments Généralisée[START_REF] Hall | Generalized method of moments[END_REF]. Cependant, le problème d'optimisation à résoudre est non-convexe, le résultat est donc une solution approchée au problème exact. La seconde approche est basée sur la complétion de la matrice de causalité de Hawkes à l'aide des premiers et second cumultants intégrés. La relaxation de ce problème s'écrit comme un problème d'optimisation convexe, ce qui nous permet donc d'obtenir une solution exacte au problème approché.Finalement, dans la partie III, nous appliquons la première méthode développée dans la partie II à des données de transactions haute-fréquence issues du carnet d'ordre du marché à terme Eurex aAEn de répondre à la question 6. La méthode est utilisée pour estimer les paramètres d'un processus de Hawkes 12-dimensionnel modélisant un actif et pour comprendre l'inØuence que les di érents évènements peuvent avoir les uns sur les autres. Ce modèle de carnet d'ordre est une extension naturelle du modèle 8-dimensionnel étudié dans[START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]. Nous augmentons ensuite la dimension du problème pour prendre en compte les évènements de deux actifs simultanément et discutons la dynamique jointe de ces deux actifs. Les méthodes non-paramétriques usuelles[START_REF] Bacry | Second order statistics characterization of hawkes processes and non-parametric estimation[END_REF] [RBRGTM14] cherchent à estimer les noyaux, ce qui restreint la dimension du modèle de carnet d'ordre pour des raisons de complexité. Notre méthode non-paramétrique estime seulement l'intégrale des noyaux, nécessite des calculs moins coûteux and passe mieux à l'échelle pour un nombre de noeuds plus grand ou un nomdre d'évènements plus important. Nous montrons aussi que la matrice de causalité de Hawkes fournit un résumé très riche des intéractions au sein du système, et peut donc devenir un outil puissant pour comprendre la structure sous-jacente d'un système présentant plusieurs types d'évènements.Nous avons maintenant assez d'éléments pour résumer les résultats principaux de cette thèse.De nombreux problèmes d'apprentissage statistique supervisé s'écrivent comme la minimisation d'une perte moyenne sur une distribution de données. D'après le principe de la minimisation du risque empirique, la perte moyenne est approchée par une moyenne des pertes sur les données observées, et un succès majeur a été de pouvoir exploiter la structure de somme pour concevoir des algorithmes stochastiques e caces[START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]. De tels algorithmes stochastiques permettent une extraction très e cace de la valeur des données massives. L'application de cette méthode aux données de survie à grande échelle, qu'il s'agisse de biostatistiques ou d'économie, est évidemment d'une grande importance. 1. Part I: Modèle de Cox à grande échelle Dans le chapitre I, nous passons en revue les progrès récents de l'optimisation convexe avec les algorithmes SGD (Stochastic Gradient Descent), du travail pionnier de [RM51] aux variantes récentes avec réduction de variance [DBLJ14] [XZ14] [SSZ13] [RSB12]. Nous introduisons ensuite la notion de processus ponctuel [DVJ07] qui fournit des outils clés pour la modélisation des événements i.e. horodatés et/ou des données de localisation. Nous introduisons enAEn le modèle à risques proportionnels de Cox [Dav72] qui relie la durée qui précède la réalisation d'un événement à une ou plusieurs covariables via la notion de taux de risque. Dans le chapitre II, nous présentons notre nouvel algorithme d'optimisation pour aider à ajuster le modèle de Cox à grande échelle.

	Plan
	Question 4. Comment adapter l'algorithme d'estimation de la régression de Cox lorsque le nombre de patients devient très grand ? passent à l1.1
	Question 5. Est-il possible de mesurer la causalité au sens de Hawkes sans faire d'hypothèse sur
	les fonctions noyaux ?
	AAEn de répondre positivement à la deuxième question, nous avons développé deux nouvelles
	méthodes d'estimation non paramétriques pour la causalité de Hawkes, plus rapides et qui

Motivations

La quantité de données collectées et stockées de façon électronique est très importante, et ne cesse de croître. L'utilisation d'outils d'analyse prédictive pour extraire de la valeur de ces données, qui est le coeur de ce que l'on appelle la data révolution, a fait ses preuves en astronomie

[START_REF] Feigelson | Big data in astronomy[END_REF]

, dans le e-commerce [MB + 12], pour les moteurs de recherche

[START_REF] Chen | Business intelligence and analytics: From big data to big impact[END_REF]

et bien d'autres domaines. Les institutions de santé se basent aujourd'hui aussi sur l'utilisation de données pour créer des modèles de traitement personnalisé grâce aux outils de l'analyse de survie

[START_REF] Murdoch | The inevitable application of big data to health care[END_REF]

. Une part importante de la recherche médicale cherche à comprendre les relations entre les covariables d'un patient et la durée avant l'occurence d'un évènement appelé défaillance (souvent la mort ou l'apparition d'une maladie).

Question 6. La méthode d'esimation de la causalité de Hawkes, introduite précédemment, peut-elle nous permettre d'avoir une compréhension plus précise de la dynamique d'un carnet d'ordre ?

Plan

Chacune des questions posée ci-dessus correspond à une partie de cette thèse.

Contexte sur les algorithmes SGD, les processus ponctuels et le modèle des risques proportionnels de Cox Dans

  ce chapitre, nous passons en revue les résultats classiques derrière les algorithmes de descente du gradient stochastique et ses adaptations à variance réduite. Nous introduisons ensuite le modèle des risques proportionnels de Cox.

	1.1.

1 Algorithmes de Descente de Gradient Stochastique Algorithmes SGD pour une distribution quelconque

De nombreux problèmes d'estimation dans le cadre de l'apprentissage statistique s'écrivent

  La comparaison des taux de convergence est toutefois di érente. Soit f deux fois di érentiable sur R d , µ-fortement convexe, ce qui signiAEe que les valeurs propres de la matrice hessienne r 2 f (µ) sont supérieures à µ > 0 pour tout µ 2 R d , et L-lisse, ce qui signiAEe que les valeurs propres sont inférieures à L > 0. Les taux de convergence avec d'autres hypothèses sur la fonction f se trouvent dans [B + 15]. On note µ § son minimiseur et on déAEnit le nombre de condition comme ∑ = L/µ. Le taux de convergence est déAEni pour les méthodes itératives comme une limite supérieure serrée d'une erreur prédéAEnie et est considéré comme la vitesse à laquelle l'algorithme converge. En notant µ t l'itération après les étapes t d'un algorithme itératif et considérant la di érence E f (µ t ) °f (µ § ) comme erreur, le taux de convergence de la Descente de Gradient est O(e °t /∑), tandis que celui de la Descente de Gradient Stochastique est O(∑/t ). Un taux de convergence de la forme O(e °AEt

  Chaque réalisation d'un processus de points ª peut être écrit comme ª = T ). Il peut être représenté de manière équivalente par un processus de comptage N ∑t } . La caractérisation habituelle du processus de point temporel se fait par la fonction intensité conditionnelle, qui est déAEnie comme la vitesse inAEnitésimale à laquelle les événements sont censés se produire après t , étant donné l'historique de N AEltration du processus qui code les informations disponibles jusqu'au temps t . Le processus ponctuel temporel le plus simple est le processus Poisson qui suppose que les événements arrivent à un taux constant, ce qui correspond à une fonction d'intensité constante

	mesure Dirac, n est une variable aléatoire à valeur entière et t	P i 's sont des éléments aléatoires n t i où ± est la i =1 ±
	de [0, t = R t 0 ª(s)d s = P n i =1 1 {t
	P(N	t +h °Nt = 1|F	t )	,
	h!0	h
	où F

i s avant t :

∏(t |F t ) = lim t est la

  Une autre façon de construire les processus de Hawkes est de considérer la représentation branchante suivante, voir [HO74] : les individus de type i , 1 ∑ i ∑ d , arrivent comme un processus de Poisson d'intensité µ i . Chaque individu peut avoir des enfants de tous types et la loi des enfants de type i d'un individu de type j qui est né ou a migré en t est un processus de Poisson inhomogène d'intensité ¡ i j (• °t ).Cette construction permet en outre de déAEnir et de mesurer la causalité entre noeuds d'un modèle de Hawkes, où les intégrales

	Introduction
	pondèrent les relations entre individus. Plus précisément, l'introduction du processus de
	comptage N	i √ j
		t
		g i j =	0 Z +1	¡ i j

°s)d N j s . (u) du ∏ 0 pour 1 ∑ i , j ∑ d .

  fois que nous observons le processus N t pour t 2 [0, T ], nous calculons les cumulants intégrés empiriques sur les fenêtres[°H ], et minimisons la di érence quadratique L T entre les cumulants théoriques et les cumulants empiriques. Nous avons prouvé la consistence de notre estimateur dans la limite T ! 1, une fois que la séquence (H T ) satisfait à certaines conditions. Notre problème peut être considéré comme une méthode généralisée des moments[START_REF] Hall | Generalized method of moments[END_REF].Pour prouver la consistence des cumulants intégrés empiriques, nous avons besoin de l'hypothèse suivante :Assumption 2. Le moitié de la taille du support du domaine d'intégration satisfait H La partie numérique, à la fois sur des ensembles de données simulées et réelles, donne des résultats très satisfaisants. Nous avons d'abord simulé des données d'événements, en utilisant l'algorithme de thinning de[START_REF] Ogata | On lewis' simulation method for point processes[END_REF], avec des formes de noyaux très di érentesexponentielle, loi de puissance et rectangulaire -et récupérons la valeur réelle du symbole G pour chaque type de noyau. Notre méthode est, à notre connaissance, la plus robuste en ce qui concerne la forme des noyaux. Nous avons ensuite appliqué notre méthode sur les 100 sites Web les plus cités de la base de données MemeTracker et sur les données du carnet d'ordres AEnanciers : nous avons surpassé les méthodes de pointe appliquée à MemeTracker et nous avons extrait des caractéristiques intéressantes et interprétables des données AEnancières. Mentionnons également que notre méthode est signiAEcativement plus rapide (environ 50 fois plus rapide) puisque les méthodes précédentes visaient à estimer des fonctions alors que nous nous concentrons uniquement sur leurs intégrales.La simplicité de la méthode, qui associe une liste de temps à une carte de causalité entre les noeuds, et sa cohérence statistique, nous a incité à concevoir de nouveaux modèles Introduction de processus ponctuels de carnet d'ordres et pour mieux comprendre sa dynamique. Les caractéristiques extraites à l'aide de notre méthode ont une interprétation économique très naturelle. C'est le but principal de la Partie III. 'approche précédente basée sur la méthode généralisée des moments a besoin des trois premiers cumulants pour obtenir su samment d'informations à partir des données pour récupérer les entrées d 2 de G. En supposant que la matrice G a une certaine structure, nous pouvons nous débarrasser du cumulant du troisième ordre et concevoir une autre méthode d'estimation en utilisant seulement les deux premiers cumulants intégrés. De plus, le problème d'optimisation qui en résulte est convexe, au contraire de la minimisation de L T ci-dessus, ce qui permet la convergence vers le minimum global. La matrice que nous voulons estimer minimise un critère simple f convexe, typiquement une norme, tout en étant cohérent avec les deux premiers cumulants intégrés empiriques.Notre problème se formule comme un problème d'optimisation sous contraintes :

	2.3 Approche par optimisation sous contraintes
	∂ °1 P °°°°! T !1	G

T , H T T ! 1 et

H 2 T /T ! 0.

Nous prouvons dans le chapitre III le théorème de consistence suivant : Result 2. Sous les hypothèses 1 et 2, la séquence d'estimateurs déAEnis par la minimisation de L T (R) converge en probabilité vers la vraie valeur G: b G T = I °µarg min R2£ L T (R) L

  T °t , L + t , L °t ,C + t ,C °t , T a t , T b t , L a t , L b t ,C a t ,C b t )où chaque dimension compte le nombre d'évènements antérieurs au temps t :

  où les dimensions P + et P °comptent les mouvements du mid-price à la hausse (baisse) due à ordre quelconque.Nous avons comparé deux couples d'actifs qui partagent les mêmes facteurs de risque. Le principal résultat empirique de cette étude concerne le couple (DAX, EURO STOXX) pour lequel les variations de prix et les variations de liquidité sur le DAX (petite tick) inØuencent principalement la liquidité sur l'EURO STOXX (grande tique), tandis que les variations de prix et les variations de liquidité sur l'EURO STOXX tendent à déclencher des mouvements de prix sur le DAX. Nous avons exécuté la procédure d'estimation sur le modèle 16-dimensionnel, nous concentrons notre discussion sur les deux sous-matrices non diagonales 8 £ 8 sur la Figure A.2, qui correspondent à l'interaction entre les actifs -l'indice D représente DAX et X pour EURO STOXX.

	a t , T b t , L a t , L b t ,C a t ,C b t )

t = (P + t , P °t , T

Part III: Capture order book dynamics with Hawkes processes

There are other types of censoring. For instance, left-censoring means the patient died or left the study before being observed : neglecting left-censoring will lead to overestimation of the survival time.

https://github.com/X-DataInitiative/tick

https://www.memetracker.org/data.html

i.e. buy orders that are executed and removed from the list

i.e. buy orders added to the list

i.e. the number of times a limit order at the ask is canceled: in our dataset, almost 95% of limit orders are canceled before execution.

A proper convex function f is a convex function taking values on the extended real line such that f (x) > °1 for all x and f (x) < +1 for at least one x.

A proper convex function is closed if and only if it is lower semi-continuous.

https://github.com/X-DataInitiative/tick

http://www.quanthouse.com

Note that we use the very same dataset as in[START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]

As was done in[START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF], for the estimation of the covariance density we take a linearly spaced grid at short time lags (until a lag of 1ms) and we switch to a log-spaced one for longer time lags. This allows to estimate the covariance on several orders of magnitude in time.

Remerciements Acknowledgements

This work beneAEted from the support of the chair "Changing markets", CMAP École Polytechnique and École Polytechnique fund raising -Data Science Initiative.

The authors want to thank Marcello Rambaldi for fruitful discussions on order book data's experiments.

Acknowledgments

This research beneAEted from the support of the Chair "Changing Markets", under the aegis of Louis Bachelier Finance and Sustainable Growth laboratory, a joint initiative of Ecole Polytechnique, Université d'Evry Val d'Essonne and Fédération Bancaire Française and from the chair of the Risk Foundation: Quantitative Management Initiative.

CHAPTER III

Generalized Method of Moments approach Abstract

We design a new nonparametric method that allows one to estimate the matrix of integrated kernels of a multivariate Hawkes process. This matrix not only encodes the mutual inØuences of each node of the process, but also disentangles the causality relationships between them. Our approach is the AErst that leads to an estimation of this matrix without any parametric modeling and estimation of the kernels themselves. As a consequence, it can give an estimation of causality relationships between nodes (or users), based on their activity timestamps (on a social network for instance), without knowing or estimating the shape of the activities lifetime. For that purpose, we introduce a moment matching method that AEts the second-order and the third-order integrated cumulants of the process. A theoretical analysis allows us to prove that this new estimation technique is consistent. Moreover, we show, on numerical experiments, that our approach is indeed very robust with respect to the shape of the kernels and gives appealing results on the MemeTracker database and on AEnancial order book data.

Keywords. Hawkes Process, Causality Inference, Cumulants, Generalized Method of Moments

Part III

Capture order book dynamics with Hawkes processes

same constant values on these blocks. Three di erent Ø 0 , Ø 1 and Ø 2 are used in the di erent blocks, with Ø 2 /Ø 1 = Ø 1 /Ø 0 = 10 and Ø 0 = 0.1. The number of events is roughly equal to 10 5 on average over the nodes. We thus obtain two datasets, the AErst one referred to as Rect10 corresponding to the rectangular kernels and the second one referred to as PLaw10 corresponding to the power law kernels. We run on these two datsets the NPHC algorithm and the ADM4 algorithm from [START_REF] Zhou | Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes[END_REF], which calibrates a single exponential kernel t ! AEØe °Øt with constant Ø, and for which we provided the intermediate true value Ø = Ø 1 . The results are shown in Figure V.2. These AEgures clearly illustrate that parametric methods can lead to very poor results when the parametrization does not represent well the data, while NPHC method gives better solutions without knowing scaling parameters Ø.

Single-asset model

In this section we apply the NPHC method to high-frequency AEnancial data. First we describe our dataset, then we compare the results of the NPHC method with those obtained with the Wiener-Hopf method of [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF] on the 8-dimensional model of single asset level-I book order events proposed in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF]. We AEnally discuss the NPHC estimation of the norm matrix associated with a "complete version" (i.e. 12-dimensional) of this model.

Data

In this paper we use level-I order book data provided by QuantHouse EUROPE/ASIA 3 for four future contracts traded on the Eurex exchange, namely the futures on the DAX and Euro Stoxx 50 equity indices, and the Bund and Bobl futures. The DAX and Euro Stoxx 50 indices track the largest stock by market capitalization in Germany and the Euro area respectively, while the Bund and Bobl are German interest rate futures on the 8.5 -10.5 years and the 4.5-5.5 years horizon respectively. The data span a period of 338 trading days from July 2013 to October 2014. For each asset, a line with the current status of the AErst levels of the order book is added to the database every time there is a change (price, volume or both). Moreover, an additional line is added in the case the change is caused by a market

Revising the 8-dimensional mono-asset model of [BJM16] : A sanity check

In [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF][START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF], the authors outlined a method for non-parametric estimation of the Hawkes kernel functions based the inAEnitesimal covariance density and the numerical solution of a Wiener-Hopf system of integral equations that links the covariance matrix and the kernel matrix. Their method has been applied to high-frequency AEnancial data in [START_REF] Bacry | Hawkes model for price and trades high-frequency dynamics[END_REF][START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF], and [START_REF] Rambaldi | The role of volume in order book dynamics: a multivariate hawkes process analysis[END_REF]. The aim of this section is to compare the newly proposed NPHC methodology with the Wiener-Hopf method mentioned above in order to assess the reliability of the new NPHC method. To this end, we reproduce the results obtained in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF].

As it was done there, we consider the DAX and Bund futures data 4 and for each asset we separate Level-I order book events into 8 categories as deAEned above: P + , P °, T a , T b , L a , L b , and C a , C b . Note that here a price move can be of any type. We then consider the timestamp associated with all events as a realization of a 8-dimensional Hawkes process and we use both the NPHC method outlined in Section 3 and the Wiener-Hopf method of [START_REF] Bacry | First-and second-order statistics characterization of hawkes processes and non-parametric estimation[END_REF] to estimate the integrated kernel interaction matrix G from the data. For the Wiener-Hopf method, we follow the same procedure as [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] and in particular we estimate the covariance density up to a maximum lag of º 1000s using a log-linear spaced grid 5 , while for the NPHC method we follow the steps outlined in Section 3 and we AEx H = 500s so to be on a comparable scale with the Wiener-Hopf method. Let us note that this scale is several orders of magnitude larger than the typical inter-event time. Indeed, on the assets considered median inter-event times are of the order of 300µs (the mean being º 50ms), with minimum time distances in the tens of microseconds.

In Figure V.3, we compare the kernel integral matrices G obtained with the NPHC method (left) with those obtained with the Wiener-Hopf approach (right) on the DAX future. Although the precise values of the matrix entries di er somewhat, as it is di cult to tune the estimation parameters of the two methods as to produce the exact same numerical results, we note that the two methods produce very consistent results. Indeed, they recover the same interaction structure and thus lead to the same interpretation of the underlying system dynamics. In our view, this represents a good sanity check for the proposed NPHC methodology. Analogous results are obtained for the Bund future. Let us also point out that the small asymmetries between symmetric interactions (such as e.g. T + ! T °and T °! T + ) can be used get a rough measure of the estimation error. In the case presented here, the average absolute di erence between symmetric interactions kernels is 0.03, which means relative error of a few percent on the most relevant interactions.

We do not comment here the features emerging from the kernel norm matrices presented in this section since they have been already discussed at length in [START_REF] Bacry | Estimation of slowly decreasing hawkes kernels: application to high-frequency order book dynamics[END_REF] and some of them will be further discussed in the next sections. Instead, here we highlight that the results of this section provide a strong case for the use of the NPHC method over the Wiener-Hopf 1 Origin of the scaling coe cient ∑ Following the theory of GMM, we denote m(X , µ) a function of the data, where X is distributed with respect to a distribution P µ 0 , which satisAEes the moment conditions g (µ) = E[m(X , µ)] = 0 if and only if µ = µ 0 , the parameter µ 0 being the ground truth. For x 1 ,..., x N observed copies of X , we denote b g i (µ) = m(x i , µ), the usual choice of weighting matrix is