{'fraction_non_alphanumeric': 0.05465506755863455, 'fraction_numerical': 0.01946111742774796, 'mean_word_length': 4.031195955833444, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-frequency Estimation of the L\'evy-driven Graph Ornstein-Uhlenbeck process

We consider the Graph Ornstein-Uhlenbeck (GrOU) process observed on a non-uniform discrete time grid and introduce discretised maximum likelihood estimators with parameters specific to the whole graph or specific to each component, or node. Under a high-frequency sampling scheme, we study the asymptotic behaviour of those estimators as the mesh size of the observation grid goes to zero. We prove two stable central limit theorems to the same distribution as in the continuously-observed case under both finite and infinite jump activity for the L\'evy driving noise. When a graph structure is not explicitly available, the stable convergence allows to consider purpose-specific sparse inference procedures, i.e. pruning, on the edges themselves in parallel to the GrOU inference and preserve its asymptotic properties. We apply the new estimators to wind capacity factor measurements, i.e. the ratio between the wind power produced locally compared to its rated peak power, across fifty locations in Northern Spain and Portugal. We show the superiority of those estimators compared to the standard least squares estimator through a simulation study extending known univariate results across graph configurations, noise types and amplitudes.


Introduction
Ornstein-Uhlenbeck (OU) models, and in particular those driven by Lévy processes, form a class of continuous-time processes with a broad range of applications, e.g. in finance for pairs trading [29,24] and volatility modelling [6,44], in electricity management [34] or neuroscience [42]. On the other hand, the availability of high-dimensional time series datasets gave rise to sparse inference for OU-type processes [12,26,40] as a way to control interactions within complex systems.
In this article, we consider the multivariate Lévy-driven Graph Ornstein-Uhlenbeck (GrOU) process on a given graph structure and the maximum likelihood estimators (MLE) introduced in [21]. This process was originally designed to quantify OU-type relationships between nodes of a continuously-observed graph. However, in this article, we suppose that observations are placed on a discrete time grid as an extension of [52] and we derive the asymptotic distributions of those discretised MLEs under a high-frequency sampling scheme [36]. Using a jump-filtering approach, we also discretise the continuous-time results from [21] to prove the stable convergence (in distribution) of the estimators to the same Gaussian distribution as in the continuous-time case when the mesh of the observation grid shrinks to zero. The stable convergence allows to use a complementary inference procedure to obtain a sparse structure, if not directly available (as in [40]), in parallel to the GrOU inference.
OU processes observed on a discrete uniformly-spaced time grid are multivariate AR(1) processes and have been well studied in this case: a least squares (LS) estimator [25,30] as well as other inference strategies such as momentbased estimators [49] or specific to nonnegative noise increments [15] have been considered. Similarly to [36], we consider the case of a non-uniform time grid with a double asymptotic assumption on the sample size, the mesh size and the time horizon. This extends the results to event-based datasets [47]. As suggested by empirical analysis [3], we ought to consider noise types that allow, in short, for both finitely and infinitely many jumps in finite time intervals-called finite and infinite jump activities respectively-and we prove that asymptotic distribution result holds in both cases. As a result, the stable asymptotic convergence holds for a large collection of Lévy driving noises, may it be the sum of Wiener and compound Poisson processes for the former, or say, generalised hyperbolic Lévy motions [23] for the latter. By property of the stable convergence, this also opens the door to the stochastic and/or time-dependent modelling of the graph structure.
To broaden the applicability of the GrOU process on high-dimensional but non-graphical datasets, we extend the regularisation approach from [26,21] to the high-frequency framework. We show the asymptotic normality and consistency in variable selection of their Adaptive Lasso scheme, a sparse inference approach that infers a lean graph structure on which the GrOU model can be applied.
We also compare quantitively the discretised MLEs with the LS estimator on a network-based electricity dataset, namely, the RE-Europe dataset [32]. A simulation study shows that the MLEs are stable across noise types, graph configurations as well as superior to the LS estimator across a variety of graph configurations and driving noises of reasonable amplitudes which corroborates and extends the results from [36].
The high-frequency stable convergence is proved with a driftless Lévy noise and under second (resp. fourth) moment integrability for finite (resp. infinite) jump activity. The infinite case also requires another stability assumption: the second moment of the Lévy measure restricted to the hypercube [−v, v] d should be O(v 2−α ) for some small v > 0 where α is the Blumenthal-Getoor index [9] of the driving Lévy noise. We validate the inference power of the MLEs on a subset of the RE-Europe data set [32] where hourly electricity-related data has been collected in 2012-2014 and projected on close to 1,500 locations of the electricity network across Europe. More precisely, we focus on the hourly wind capacity factor measurements at 50 of those locations where 24 are located in Portugal and 26 are in Northern Spain along the Atlantic coast. We recall that the wind capacity factor is the proportion of energy produced compared to the peak power production at this location hence between 0 and 1. The data represents the local supply of wind energy that ought to be distributed across the networks: the GrOU model captures the geographical relationships between nodes while maintaining the graphical structure of the electricity grid on which the energy transfer take place. We will infer the model parameters and noise distributions on this dataset and use those parameters in the simulation study. This study focuses on testing the reliability of the estimation across noisy types, graph structures and mesh sizes by recovering the model parameters obtained from parametrically bootstrapped time series.
In Section 2, we recall the definition of the GrOU process when observed on a discrete time grid and introduce the corresponding discretised MLEs. In Section 3, we describe the high-frequency sampling scheme which allows the asymptotic stable convergence of the discretised GrOU MLEs given thereinwith separate treatments of the finite and infinite jump activity cases. Section 4 is devoted the definition of the Adaptive Lasso scheme and the study of its asymptotic properties. Finally, in Section 5, we apply the discretised MLEs on the subset of the RE-Europe data set [32]. Additionally, we perform the recovery of Lévy increments [16] to infer the noise distribution via likelihood maximisation or the generalised method of moments. Those fitted distributions are used in the simulation study presented in Section 6. Finally, the technical lemmas and propositions are stated in Appendix D whilst the proofs are relegated to Appendix E & F for a clearer layout of the argument. Finally, implementation details, R code and reproducible plots are available on GitHub (numerical routines 1 , end-to-end scripts 2 ).

The Graph Ornstein-Uhlenbeck process and estimators
In this section, we recall the continuous-time framework of the Graph Ornstein-Uhlenbeck (GrOU) process with two parametrisations as defined in Sections 2-4, [21].

Notations
We consider a filtered probability space (Ω, F, (F t , t ≥ 0), P 0 ) to which all stochastic processes are adapted. We consider a Lévy process (L t , t ≥ 0) (i.e. a stochastic process with stationary and independent increments and continuous in probability and L 0 = 0 d , P 0 − a.s.) [14,Remark 1] without loss of generality assumed to be càdlàg. For a stochastic process (Y t , t ≥ 0), we write Y t− := lim s↑t Y s for any t ≥ 0. For any probability measure P, we denote by P t its restriction to the σ-field F t for any t ≥ 0. Regarding the convergence of random variables, we write We denote by λ leb the one-dimensional Lebesgue measure. For a non-empty topological space, B(S) is the Borel σ-algebra on S and π is some probability measure on (S, B(S)). The collection of all Borel sets in S ×R with finite π⊗λ lebmeasure is written as B b (S × R). Also, the norms of vectors and matrices are denoted by · . We usually take the Euclidean (or Frobenius) norm but due to the equivalence between norms, our results are not norm-specific and are valid under any norm in R d or M d (R). In addition, for an invertible matrix M ∈ M d (R), we define x, y M := x M −1 y for x, y ∈ R d . Finally, for a process X t = (X (1) t , . . . , X (d) t ) ∈ R d we denote by [X] t the matrix ([X (i) , X (j) ] t ) of quadratic co-variations on [0, t] for t ≥ 0.
In this article, ⊗ denotes the Kronecker matrix product, is for the Hadamard (element-wise) matrix product and vec is the vectorisation transformation where columns are stacked on one another. We denote by vec −1 (x) := (vec(I d×d ) ⊗ I d×d )(I d×d ⊗ x) ∈ M d (R) for x ∈ R d 2 the inverse vectorisation transformation.
Then, under Assumption 1, there exists a unique strictly stationary solution to Equation (2.2.1) given by Similarly to [25], we consider a Lévy-driven OU process Y observed on a discrete time grid 0 = t 1,N < · · · < t N,N := T N where N ∈ N which leads to multivariate AR(1) process representation given by Notation 2.1. When there is no ambiguity, we use t k for t k,N and Y k for Y t k .
We study the asymptotic properties of a couple of discretised estimators for Q under two parametrisations specific to graphs which we describe in the following section.

A graphical interpretation of the OU process
The components of Y are interpreted as the nodes of a graph structure and the undirected links, or edges, between them are given by a so-called adjacency matrix (or graph topology matrix) A = (a ij ) ∈ M d ({0, 1}) where a ii = 0 for i ∈ {1, . . . , d} and, for i = j, a ij = 1 if nodes i and j are linked and 0 otherwise. We define the degree of any node i as n i := 1 ∨ j =i a ij and we also define the row-normalised adjacency matrix Assumption 2. We assume that A is deterministic, constant in time and known. Following [21], the dynamics matrix Q is parametrised in two different ways specific to graphs, i.e. they include the adjacency matrix A and model interactions within the graph.

The two-parameter formulation
The θ-GrOU parameter θ = (θ 1 , θ 2 ) ∈ R 2 give way to the parametrisation This yields single network and momentum effects for the whole graph respectively parametrised by θ 1 and θ 2 . The resulting autoregressive form is given by

A more detailed formulation
We also consider the ψ-GrOU parameter ψ ∈ R d 2 which in turn yields the parametrisation where is the Hadamard product. In this case, momentum and network effects are specific to each node and its neighbours. Similarly, the autoregressive form is given by for any i ∈ {1, . . . , d}. For this process to be well-defined, Q(ψ) ought to be positive definite [38] and we consider the following sufficient condition [

The Graph Ornstein-Uhlenbeck process
Here, we recall the definition of the Graph Ornstein-Ulhenbeck (GrOU) process in its two different configurations. Definition 2.2. The Graph Ornstein-Uhlenbeck (GrOU) process is a càdlàg process (Y t , t ≥ 0) satisfying Equation (2.2.1) with a Lévy process (L t , t ≥ 0) where Q is given by either Equation (2.3.1) or by Equation (2.3.2) such that Q is positive definite. This process is then called a θ-GrOU process or a ψ-GrOU process, respectively.
Since we consider the inference of GrOU processes observed on a discrete time grid, we first introduce the corresponding continuous-time estimators.

Continuous-time estimators
To estimate the parameters θ and ψ, [21] propose a likelihood framework under Assumptions 1 & 2 and study the corresponding maximum likelihood estimators (MLE). More precisely, for a θ-GrOU process which is observed in continuous time, has finite second moments and satisfies Assumption 3, the MLE for θ = (θ 1 , θ 2 ) ∈ R 2 is given by for any t ≥ 0, where we define with Y c , the continuous P 0 -martingale part of Y, defined by dY c t = −QY t− dt + dW t for a d-dimensional Brownian motion with covariance matrix Σ. Similarly, for a ψ-GrOU process which is observed in continuous time, has finite second moments and satisfies Assumption 4, the MLE for ψ ∈ R d 2 is given by

Discretised estimators
Recall that we denote by Y c the P 0 -martingale part of Y. Note that the parameters sets Θ and Ψ are equivalent for the θ-GrOU process and Assumption 4 boils down to Assumption 3 in this case. Therefore, in general, we only mention Assumption 4. We use the jump-filtering approach from [36] and discretise the continuous-time estimators defined in [21] Also, for a collection (v (i) N , i ∈ {1, . . . , d}) ⊂ (0, ∞), we introduce the R d 2 -valued random vectors given by: where i, j ∈ {1, . . . , d}. We define the discretised unfiltered estimator as and the discretised jump-filtered estimator by Remark 2.4. Both θ t and ψ t feature the diffusion matrix Σ ∈ S ++ d which is generally not available and we have reformulated them without Σ as in the proof of Proposition 3.3.2, [21]. Remark that as given in the proof of Proposition 3.3.2, [21]. We can rewrite such that the estimator ψ t is given by We proceed similarly for θ t .
Remark 2.5. The collection (v (i) N , i ∈ {1, . . . , d}) consists of jump-filtering threshold values: it relies on the assumption that when noise increments are small, they are most likely due to the continuous martingale part of the Lévy process (i.e. the Brownian motion part). On the other hand, larger values would correspond to the contribution of the jump part J (see Appendix A) and should be discarded as such. However, as explored in Section 3 and beyond, the exact values chosen for these thresholds are crucial as they balance the Brownian increments against the jumps. This simple methodology is instrumental for those estimators to be consistent and asymptotically normal as well as usable in practice (Sections 4 & 5).
Since the continuous martingale part is not available in practice, the discretised jump-filtered estimator ψ N is the only one usable for applications. Informally, the main result of this article is that the discretised jump-filtered estimator satisfies the stable central limit theorem given by under a collection of assumptions. We also obtain a similar estimator for θ.
Definition 2.6. Assume that Assumptions 1-3 hold for a θ-GrOU process Y. We define the θ-GrOU discretised jump-filtered estimator This estimator is directly obtained from the formulation of ψ t above since vec(Q(θ)) = θ 2 vec(I d×d ) + θ 1 vec(A). We have Note that if each node has at least one neighbour, we have ρ(A) = I { X =0} (A). We deduce that it suffices to prove the asymptotic convergence properties for the ψ-GrOU and deduce the result for the θ-GrOU process using Definition 2.6.
To convert the results involving the continuous martingale part of Y to the practical jump-filtered estimator, we split the convergence study into two sections for both finite and infinite jump activities.

Asymptotic theory with discretised high-frequency observations
We extend the high-frequency sampling scheme from Section 1, [36], to discretise the observations of the GrOU process and study the asymptotic distributions of ψ N and θ N .

High-frequency framework
The framework consists of two long-span asymptotic assumptions: Assumption 5. Let (Y t , t ≥ 0) be a GrOU process. Suppose that we are given observations Y t 1,N , . . . , Y t N,N at times 0 = t 1,N < · · · < t N,N := T N . Also, suppose that, as N → ∞, we have T N → ∞, a smaller mesh size ∆ N := max i |t i+1,N − t i,N | → 0 and a relative growth such that The framework described in Assumption 5 is said to be a high-frequency sampling scheme because ∆ N goes faster to zero than T 2 N goes to ∞. The number of samples N is curbed since N ∆ N must grow as fast as T N as N → ∞. Note that this assumption implies that data is synchronously available across all nodes or components of Y. The asynchronous case is not relevant for our application of interest, but is important for other applications; see [18,Chapters 3 & 4] and [2] for respectively autoregressive and log-Gaussian statistical approaches whilst [8] proposes a recurrent neural network approach. All of these models are applied on asynchronous financial market closing prices. However, this is beyond the scope of this article and we leave this open problem for future research.
By extending Assumption 3.1, [36] to the multivariate case, we provide the following additional sets of assumptions to obtain the asymptotic consistency and normality of the discretised estimators. First, regarding the GrOU process and the Lévy process L: Assumption 6. Let (Y t , t ≥ 0) be a GrOU process. Assume that: (i) the square integrability of the GrOU process Y, i.e. E Y t 2 < ∞ (ii) the underlying Lévy noise drift w.r.t the truncation function τ (z) := We impose an additional assumption on the high-frequency sampling scheme: Assumption 7. We assume that there exists β (i) ∈ (0, 1/2) such that

Asymptotic convergence of discretised estimators
A key result about the asymptotic distribution of the discretised unfiltered estimator ψ N is given as follows: Proof. See Appendix B.
Proof. Lemma 3.1 yields the stable convergence of T 1/2 N (ψ N − ψ) towards of Gaussian distribution with a finite and constant covariance matrix. Thus, ψ N − ψ has a O(T −1/2 N ) vanishing covariance matrix which becomes the Dirac delta distribution at zero as N → ∞. From Corollary 3.6, [28], this ensures the convergence in probability to zero of ψ N − ψ.

Jump thresholds
The jump thresholds play an important role in controlling how much of the noise increments are attributed to the continuous martingale part or as jumps (Remark 2.5). We use the thresholds defined by v (j) N , for β (j) ∈ (0, 1/2) and j ∈ {1 . . . , d}, for which the parameter β (j) allows for larger increments (as β (j) → 0) or narrows down to the usual Gaussian scale (as β (j) → 1/2) as a power of the mesh size. In a sense, the asymptotic behaviour relies on ∆ N in the theoretical study (β (j) being given as a parameter) and on β (j) in practice as it needs to be estimated (∆ N being imposed by the data at hand). The data-based aspects of β (j) and the mesh size ∆ N are discussed in Sections 5.4 & 6.4 respectively. We generalise the central limit theorem above to jump-filtered discretised estimators in both the finite and infinite jump activity cases. The technical lemmas are presented in Appendix D and proved in Appendix E.

Finite jump activity case
In this section, we consider the finite jump activity case, i.e. R d ν(dz) < ∞. The finite jump activity case is the case where the Lévy noise corresponds to the sum of a Brownian Motion and a compound Poisson process, see Appendix A. An additional assumption regarding the distribution of smaller jumps is required and reads as follows: In addition, suppose that the marginal distributions of the jump heights are such that where F (i) is the i-th node jump marginal cumulative distribution function for i ∈ {1, . . . , d}, To obtain confidence intervals on the GrOU parameter estimation, we prove a central limit theorem for the discretised jump-filtered ψ-GrOU estimator under finite jump activity as follows: . . , d}, the discretised jump-filtered estimator ψ N satisfies Proof. By Slutsky's lemma and Lemmas 3.1 & 3.6.
Next, from its definition, the θ-GrOU parameters follows a similar central limit theorem.
. . , d}, the discretised jump-filtered estimator ψ N satisfies where G ∞ is defined as Proof. This follows from Theorem 3.3 and the definition of θ N .
Remark 3.5. Under Assumptions 6 and 8, the Lévy-Itô decomposition (see Appendix A) yields that there exists a centred Gaussian Lévy process W t with covariance matrix Σ and almost-surely continuous paths as well as a compound Poisson process J t independent of W t such that Indeed, there exists a family of Poisson processes N (j) t with elementwise intensity To prove those afore-mentioned results, we rely on the consistency of the jump-filtered estimator ψ N with respect to the unfiltered estimator ψ N given in the following lemma: Proof. See Appendix E.1.
From Corollary 3.2, this means that ψ N is also consistent.

Infinite jump activity case
The infinite activity models have been investigated numerous times as they provide realistic scenarios in financial applications, avoiding for instance to introduce a Brownian component as empirically determined in [17]. More precisely, it entails assuming that R d ν(dz) = ∞. The Lévy-Itô decomposition with respect to τ (z) := I {x∈R d : x ≤1} (z) describes the pure-jump process J as the sum of jumps larger than one J 1 and jumps smaller than one J 2 , such that J d = J 1 + J 2 . Moreover, the infinite jump activity case features an infinite number of those small jumps in any time interval. In this context, an additional compensation on the Lévy measure ν is provided to control J 2 , see Appendix A. We also extend Assumptions 4.1, [36] to our framework as follows: Assumption 9. For a GrOU process Y, assume that as well as the following properties: (i) Y has bounded fourth-order moments, i.e. for any t ≥ 0, E Y t 4 < ∞. (ii) The squared magnitude of small jumps remains controlled in the sense that there exists α ∈ (0, 2) such that, there exist 0 > 0 and n 0 ∈ N such that for all ≤ 0 and n ≥ n 0 , one has for all k ∈ {0, . . . , N − 1} and i ∈ {1, . . . , d}.
The relationship between the small jumps behaviour between the finite and infinite cases is given by Remark 4.3, [36]: if the Lévy measure ν is finite, then Assumption 8 implies Assumption 9-(ii) if , for all i ∈ {1, . . . , d}. This means that a small-jump density that is regular enough around zero is sufficient to imply the jump-related regularity conditions. In addition, Assumption 9-(iii) yields that small jumps described by the Lévy measure of J 2 (i.e. ν restricted to [−1, 1]) are centred around zero. A related result is that if the Lévy measure of J 2,(j) is symmetric (potentially restricted to an interval (− 0 , 0 )) then Assumption 9-(iii) holds [36,Lemma 4.4]. This includes symmetric αstable distributions or zero-mean generalised hyperbolic distributions [23] such as the Student's t-distribution, the Variance-Gamma distribution or zero-mean Normal-inverse Gaussian distribution (see Section 5).
We hereby propose to adapt the proofs of Section 4, [36] to the multivariate framework. As a stepping stone for the consistency and asymptotic normality of the jump-filtered discretised estimators, we prove the consistency of the jumpfiltered estimator with respect to the unfiltered one similar to Lemma 3.6 as follows: Proof. See Appendix F.9.
. . , d}, the discretised jump-filtered estimator ψ N satisfies Proof. By Lemma 3.1, we know that the discretised unfiltered estimator ψ N converges stably in distribution as follows: By Lemma 3.7 we have that T We conclude by Slutsky's lemma.
Similarly, for the θ-GrOU estimator, we obtain the following theorem . . , d}, the θ-GrOU discretised jump-filtered estimator θ N satisfies Proof. By Theorem 3.8 and the definition of θ N .
As noted above in the finite activity case, from Corollary 3.2, this means that ψ N is also consistent.
We have now defined the discretised estimators in Section 2 and provided central limit theorems for both the θ-GrOU and ψ-GrOU estimators in Section 3 for Lévy driving noises for both finite and infinite jump activities. In the next section, we describe the different methodologies used to infer the GrOU parameters as well as recover the Lévy increments and estimate their distribution as means to perform parametric bootstrap.

Asymptotic convergence of the Adaptive Lasso regularisation
Both the θ-GrOU and ψ-GrOU demand the graph adjacency matrix A to be fully specified (Eqs. (2.3.1) and (2.3.2)), which restricts the applicability of those processes to data that is based solely on graphs. Sparse inference has risen as a way to control the dimensionality of the parameter vector by providing lean graph-like structure in high-dimensional problems [19,35].
In this section, we extend the L 1 -penalised likelihood function from [26,21] to discretely-observed data to estimate a candidate adjacency matrix: (i) which is fully determined by the information given in the dataset; (ii) whose sparsity can be tuned through an explicit parameter. This Lasso-type regularisation scheme is shown to estimate the adjacency matrix consistently and we provide stable asymptotic normality of the dynamics matrix estimator.

Adaptive Lasso
For a Lévy-driven OU process which satisfies Eq. (2.2.1) with a true but unknown dynamics matrix Q 0 ∈ S ++ d , a Lasso regularisation which penalises the likelihood with the L 1 -norm of the dynamics matrix cannot achieve both aforementioned properties in, say, the Gaussian case [53]. An Adaptive Lasso (AL) procedure provides superior theoretical properties as shown in [26] with a Brownian motion as the driving noise and in [21] with a driving Lévy process with continuous-time observations. We solely extend the latter to the high-frequency framework of Section 3.1. For fixed parameters λ ≥ 0 and γ > 0, the AL scheme is defined as where denotes the Hadamard product and where the absolute value in the denominator of the penalty is evaluated elementwise. The discretised log-likelihood N (·) (equivalent to Eq. 6, [21], which is in continuous time) is given by which penalises more those that are expected to be zero. We also recall that for i, j ∈ {1, . . . , d}, we have:

Asymptotic properties
We prove that the AL scheme is consistent in variable selection as well as asymptotically normal under the existence of an oracle that knows the true dynamics matrix Q 0 . Note that the penalty parameter λ is implicitly a function of the time horizon T N which we write λ = λ(T N ). The support of vector or matrix is defined as follows: Notation 4.1. The support of a vector or a matrix x is denoted supp(x) and is defined as the set of indices of non-null coordinates of x. For a given set of indices I, we denote by x |I , the restriction of x to the indices in I and x |I×I to the indices in I × I.

Application to wind capacity factor measurements
In this section, we present the RE-Europe data set on which we apply the discretised GrOU process whose inference serves as a foundation for the simulation study.

Data set and motivation
The RE-Europe dataset [32] offers an approximation of the hourly wind capacity factor (between 0 and 1) across 1,500 locations using the COSMO-REA6 data set [11] for the years 2012-2014. This data set consists of a grid of wind speed measurements linearly interpolated to obtain wind speed at turbine height (80 m) with a 7 × 7 kms 2 spatial resolution. It is then further processed to obtain wind capacity factors and projected on the 1,494 locations of the European electricity network [32, page 7] of which we consider a small portion. Wind capacity factors represent the proportion of the maximum power generation that is generated at a given time. We focus on a subset of 50 locations where 24 are located in Portugal and 26 are in Northern Spain along the Atlantic coast (see LHS, Figure 1). This 50-node network reflects how distributed power generation capabilities would behave in reality with operational constraints: four subnetworks are disconnected from one another and power demand needs to be fulfilled by the local generation of power. The GrOU process describes the geographical dependencies between the nodes while maintaining the graphical information of neighbours and disconnected subset of nodes.
The 50-node network and the related wind capacity factors data will be referred to as the RE-Europe 50 dataset which comprises of 25,000 hourly data points.

Methodology
For the application on the wind capacity factors data set in this section and the simulation study presented respectively in Section 6, estimating the model parameters is only one of the necessary steps. We devote this section to clarifying the steps involved in preparing the data, obtaining the Lévy increments recovery and inferring the increments' distribution for parametric bootstrapping.

Algorithm
For applications and the parametric bootstrap procedure, we proceed through the following steps 1. Separate trend, seasonal and remainder time series using a LOESS regression or alternative methods [33]; 2. Fit either the θ-GrOU or ψ-GrOU model with the MLE to the remainder time series; 3. Recover the Lévy increments which we denote (∆ k L : k ∈ {1, . . . , N − 1}) (Section 5.2.3); 4. Infer the Lévy driving noise distribution from those increments (Section 5.2.4).
In the following sections, we give details on the steps 2 through 4. Note that Step 4 relies on the knowledge of adjacency matrix A.

With an unknown adjacency matrix
If A is not known, one can proceed in two different ways: (i) by inferring it on the same dataset using a particular scheme, such as Lasso regularisation (see Section 4) in lieu of Step 2 above; (ii) by estimating a sparse adjacency matrix on a potentially different dataset with a consistent estimator where A is either a deterministic or random matrix. Then, we estimate the dynamics as follows If A is deterministic then the pairwise convergences described above would hold in distribution by a standard argument [28, p. 6-7]. On the other hand, if A is a random variable, the pairwise convergence in distribution will not even hold in general without the stability [28, Example In this sense, the stable convergence in distribution paves the way towards an extension to a random and/or time-dependent adjacency matrix (e.g. as an additional, potentially non-ergodic, process) [51,5].

Lévy increments recovery
In this section, we describe the recovery of Lévy increments to infer their distribution. [16] approximate the noise increments from the data and the estimated matrix Q, i.e.
Since Q is not available in practice, we use the estimates Q( θ N ) or Q( ψ N ) and, for the integral, the k-th absolute moment of the approximation error is asymptotically O(∆ 2k N ) provided that the process is observed on a uniformlyspaced time grid with increment ∆ N and that the k-th moment of the Lévy process is finite [16,Proposition 5.4]. We adapt to non-uniformly-spaced time intervals given our assumption of an asymptotically zero mesh size ∆ N as N → ∞ and to estimate the 1-step increment, we approximate it crudely by Those estimated, or recovered, Lévy increments are denoted ∆ k L, for k ∈ {1, . . . , N − 1}.

Lévy driving noise inference
Once the increments are recovered, we proceed to fit those to a finite and an infinite jump activity distribution. For the former, we first leverage a generalised method of moments inference, as inspired by [42], to estimate the Brownian motion covariance matrix. Jumps are either modelled using a d-dimensional Gaussian distribution or a generalised hyperbolic distribution (GHYP), see below. For the latter, we directly model the noise as generalised hyperbolic Lévy motions [23].
Generalised hyperbolic distribution It is a Normal mean-variance mixture distribution which features semi-heavy tails along with an explicit likelihood function. Thus, it is particularly suitable for jump distribution modelling. We assume that the jumps can be decomposed as follows . We use an Expectation Maximisation algorithm called the Multi-Cycle Expectation Conditional Maximization (MCECM) algorithm [41], as implemented in the R package ghyp, to fit the full 6-parameter GHYP density directly on increments.
Finite jump activity case Following Remark 3.5, we assume that L is the sum of a correlated d-dimensional Wiener process W with covariance matrix Σ and a d-dimensional compound Poisson process with finite intensity λ > 0 and zero-mean Gaussian or GHYP jumps. Gaussian jumps are modelled by a d-dimensional Gaussian random vector N (0 d , Σ J ) where Σ J ∈ S + d , similarly to the GHYP notation. In the spirit of [42,27,3], we calibrate both the correlation matrix of the Wiener process and the compound Poisson process with Gaussian jumps using the realised (co)variance matrices. In this case, for i, j ∈ {1, . . . , d}, we define the (i, j)th component of the realised covariance by The jump-filtered equivalent is defined for a unique fixed β ∈ (0, 1/2) as follows for some η > 0. That is, increments vectors are considered to be jumps if the L 2 -norm of the recovered increments are above η · ∆ β N : for instance, η = 1 implies that jumps include all increment vectors such that |∆ k L (j) | > ∆ β for some j ∈ {1, . . . , d}. Remark 5.1. Note that (5.2.1) relies on the assumption of a fixed volatility (i.e. a fixed η > 0), e.g. a fixed Gaussian covariance structure. An extension of this estimator to a variable volatility (similarly to [7]) would be given by where σ k is a spot volatility estimate at time t k for k ∈ {0, . . . , N − 1}. Remark 5.2. Alternatively to estimate RV(L) C N , one could have used the bipower variation [7] Finally, the realised (co)variance of the jump part defined as Since ∆ k L P − − → ∆ k L as ∆ N → 0, those quantities approximate quadratic covariation matrices as ∆ N → 0 stably in distribution [27,4]. Hence, we have for N large enough Σ ≈ RV( L) C N /T N , whilst the jump intensity λ is estimated by likelihood maximisation given that the probability of n jumps in [0, T N ] is proportional to Recall that we assume that jumps occur when ∆ k L 2 > η · ∆ β N . Finally, for Gaussian jumps, we obtain directly Infinite jump activity case We consider the generalised hyperbolic Lévy motions (GHYP) for their important distributional flexibility to accommodate for fat tails [23] and use the likelihood maximisation directly on the recovered Lévy increments without filtering for jumps. Remark 5.3. The parametric bootstrap used in the following sections is supported by the theoretical guarantees presented in [48,1]: consistent parameter estimators for the driving Lévy noise distribution are enough to replicate the true dynamics of the GrOU process. Although those results are proved for distributions characterised by their first two moments [1, Th. 3.1], we assume it also holds for GHYP distributions.

Exploratory analysis and driving noise inference
We clean our data set using a LOESS regression (R package stl) with daily seasonality and a daily trend rolling window. The output has no unit root (via the Dicker-Fuller test with 99.9% significance on each marginal) and has zero mean.
Remark 5.4. The wind capacity factor values are restricted to [0, 1] whilst the processed time series are unrestricted real values.
Using the RE-Europe 50 data set, the θ-GrOU parameters is estimated by with (infinite activity) standard deviations in parenthesis on 500 samples of 12,500 observations each through a parametric bootstrap. Note that the finite activity standard deviations are similar up to the second decimal hence omitted. Equation (5.3.1) yields that, on average, each node's value acts negatively on its increments, i.e. they have a negative momentum, whilst the neighbouring nodes act positively on the increments as they sort of aggregate the neighbours' capacity towards the node. That is, if the capacity of a node spikes, it may indicate that neighbouring nodes might increase soon after that. This said, the momentum term in absolute value is four times as big as the network term which indicates the strong autoregressive nature of wind flows. The ψ-GrOU parameters given on the right-hand side in Figure 1 show that Portuguese locations seem to have lower momentum than Spanish ones (top violin plot) whilst the negative Spanish network effects can be larger in absolute values than the Portuguese counter parts. This said, both the momentum and networks effect means cannot be distinguished between the two countries with a two-sample t-test with a p-value larger than 90%.

Fitting the recovered increments
Noise marginals are symmetric and have fatter tails than a single Gaussian distributions (see Figure 2) hence the sum of two Gaussian marginals with a standard deviation larger for jumps or GHYP distributions are adapted to this data set.
Regarding the finite activity modelling, we first ought to specific the jump cutoff value η * such that the recovery increments vectors ∆ k L are considered jumps if ∆ k L 2 > η * ∆ β N . We discuss the value of β in Section 6.3 and use β = 0.499. In Figure 3, we present the empirical probability of having a jump for different values of η along with the variance captured by either the continuous part (through Σ) or the jump part (through λ · Σ J ) under Gaussian jumps. Covariance matrices are estimated as in Section 5.2.4 and we compute the average ratio between the diagonal elements of Σ and λ · Σ J and that of the empirical covariance matrix of (∆ k L : k ∈ {1, . . . , N }), respectively, for various values of η. We observe that we obtain the coverage equilibrium for η = 1.7, Choosing the jump cutoff value η: (left) empirical probability that ∆ k L 2 /∆ β N > η; (right) Variance coverage of the continuous part and jump parts with 95% error bars. Lines corresponding to η = 1 (50% quantile of jump norms) and η = 1.7 (80% quantile of jump norms) are drawn.
which corresponds to the 80% quantile of the L 2 -norm of the increments, and we set η * = 1.7.
We find that noise is mostly clustered by country ( Figure 4): we have higher correlations between nodes of the same country (top-left and bottom-right quadrants) and lower correlations for pairs of cross-country nodes. Correlations are positively skewed with minima of −40% and mostly positive (between +20% and +60%) or close to zero on the second and third quadrants. The different topologies between Spain (multiple subgraphs) and Portugal (a single 24-node connected graph, see Fig. 7   As a check, we generate samples from the fitted distributions and compare their quantiles to the standardised empirical quantiles in Figure 5. It features the Q-Q plots of finite activity jump distributions of a compound Poisson with Gaussian or GHYP jumps, and the infinite activity GHYP distribution. Gaussian jumps fall short in terms of tail heaviness as shown by the S-shaped Q-Q plots. The GHYP distribution perform better in that respect as shown by the other two the GHYP-based samples. We choose to only use the compound Poisson with GHYP jumps in the rest of the study. We note that across the three plots, a slight bump is observed for quantiles between 0 and 1 which may indicate that the skewness of the fitted distributions is not sufficient to represent the data correctly.

Simulation study
In this section, we assess the reliability of the inference approach on data simulated from the model itself. By avoiding the potential model misspecification, we study the impact of the graph topology, the jump threshold parameters and we compare the MLE to the least squares (LS) estimator. We use the empirical parameters values θ N found in Section 5 as true. We perform parametric bootstraps of our model on paths of N = 12, 500 observations each which corresponds to 50% of the original data set length and on a uniform time grid (k∆ N : k ∈ {0, . . . , N − 1}) where ∆ N = 1/12. For the finite activity driving noise, we generate the median number of jumps across dimensions in each dimension given we work on a fixed time interval [0, T N ] to ensure a reliable number of jumps are being drawn.

Empirical validation
We check that the θ-GrOU estimators converge the true value as the sample size grows and that they each follow a univariate Gaussian distribution by performing the inference on only the first N = 500; 4, 500; 8, 500; 12, 500 observations of each path. We restrict ourselves to the infinite jump activity case since it is the only example with a non-Gaussian driving noise; we also only work with the RE-Europe 50 network for simplicity. Gaussian Q-Q plots and histograms over the collection of 500 independent fitted values for N = 500 and N = 12, 500 are presented in Figure 6. The Q-Q plots for both θ 1 and θ 2 are straight and the Gaussian fit on the histogram has the right width. This said, with limited data (N = 500), the Gaussian mode is lower than the true value θ 1 = −1.549 and greater than the true value θ 2 = 5.525 (by a number of standard deviations for the latter). The Gaussian fit is particularly suited for θ 1 with N = 12, 500. Table 6.1 features the numerical values of the experiment: as N grows, the bias (and percent bias) reduces as expected although the fit with only 500 samples is already satisfactory with 3.73% and 11.94% biases respectively for θ 1 and θ 2 . Those biases boils down to a few tenths of a percent for N = 12, 500. The bias and standard deviations both shrink by a factor of 50 between the cases N = 500 and N = 12, 500. Their rates of convergence against the sample size N are estimated to be 1.23(.106) and 0.55(.006), respectively, using a standard linear regression on the log-log scale; note that both are significant at the 5% level. This corresponds to the expected rate of 1/2 described in Section 3.

Impact of the graph topology
In this section, we study the impact of the adjacency matrix, or graph topology, A. From Equation (2.2.1), it seems reasonable to postulate that this said topology impacts the inference performance of the estimators-especially for the off-diagonal parameter θ 1 . The data corroborates this effect which we interpret as a consequence of the value dissipation across a higher number of neighbours; formally, its expression features the quadratic impact of synchronised peaks between two variables. Figure 7 presents three artificial network configuration by increasing order of connectivity (here shown by some key statistics of the node degrees in parenthesis): Polymer (50 nodes in a straight line, mean = 2.96, median = 3, min = 1, max = 2), Lattice (a grid of 9 × 9 nodes where nodes are connected to their closest neighbours and one additional node linked to only one node of the grid, mean = 3.36, median = 4, min = 2, max = 4), Complete (a 50-node complete graph, mean = 50, median = 50, min = 50, max = 50), along with the RE-Europe 50 network configuration (mean = 1.12, median = 1, min = 0, max = 4).

Polymer
Lattice Complete RE−Europe 50 In Figure 8, violin plots of the inference on 100 paths with 12,500 samples each are represented for those four configurations.
The MLE estimator performs well across all graph configurations except for the off-diagonal parameter θ 1 for the complete graph which is substantially noisier. Quantitatively, the row-normalised parameter θ 1 /n i is low compared to θ 2 where n i = 50 which induces some numerical instability. Qualitatively, the complete graph dissipates the signal at a given node to all its neighbours whilst receiving a small portion from each of its neighbours' signal which causes some instability. The performance for θ 2 is similar across the examples since each node is equipped with a self-loop in all graph configurations.

Tuning the jump threshold exponent β
An important set of parameters are the jump thresholds (β (i) : i ∈ {1, . . . , d}) ⊂ (0, 1/2) which we assume to be equal across the marginals given they are centred and with similar standard deviations (≈ 5%). A jump threshold too large, i.e. β → 0, (resp. too small, i.e. β → 1/2) would lead to underestimate (resp. overestimate) the frequency of jumps and overestimate (resp. underestimate) their amplitudes (Remark 2.5). Parametrically bootstrapping 500 paths using the estimated infinite activity driving noise distribution led to the bias study presented in Figure 9. As β → 1/2, the bias for the off-diagonal term θ 1 decreases in all four cases with a maximum change of 0.005 in absolute value. It becomes negligible for RE-Europe 50 and Complete topologies whilst it gets even more negative for the Complete and Lattice topologies. However, the changes in bias for θ 1 are substantially smaller than the biases for θ 2 : the biases change by a higher order of magnitude from approximately −0.025 to −0.010 for all four topologies.
This contrasts with the increased variability in the estimation for θ 1 compared to θ 2 depending on the topology (Section 6.2). Therefore, we take a value of β close to 1/2 (β = 0.4999) for two reasons: (a) to limit the estimation bias for both θ 1 and θ 2 with the EU-Europe 50 topology (b) since θ 2 benefits the most from this change in β while θ 1 is already more volatile. This is in accordance with the literature [10, Section 4, p. 1749].

Comparing with the least squares estimator
In this section, we compare the high-frequency Lévy-driven OU MLE estimator to the LS estimator as studied in [25] which estimates exp(−Q∆ N ) and we extend this comparison to the multivariate context. We report the Relative where θ is given in Eq. (5.3.1) as the true parameters for the paths simulation. The smaller the REM, the better the inference strategy is at recovery the true parameters whilst allowing to compare both the MLE and the LS estimators.

Noise amplitude robustness
To measure the robustness of estimators against the noise amplitude, we multiply our fitted driving noise distribution by a noise multiplier scalar σ ∈ {0.5, 1, 5, 10, 100, 1000}: we generate our paths with the noise increments (σ∆ k L : k ∈ {1, . . . , N }). This puts the robustness of the θ-GrOU MLE to the test as the noise amplitude gets larger across graph topologies with both finite and infinity jumps activities. In Figure 10, we vary the noise multiplier parameter σ ∈ {0.5, 1, 5, 10, 100, 1000} to generate paths driven by a noise with either finite (MLE (fin)) or infinite (MLE (inf)) jump activity.
The LS estimator has a constant performance across the jump activities and graph topologies with the REM ranging from 0.37 to 0.41 in the finite and infinite cases. Thus, we only plot the finite cases. Similarly, the MLE performs equally under finite and infinite activities.
As shown [36] in the univariate case, the MLE estimator is superior to the LS estimator with the REM ranging from 0.05 to 0.20 for σ ∈ {0.5, 1, 5, 10}. As σ gets larger, the MLE performance gets poorer with the REM above 0.5 for σ ∈ {100, 1000} in both the finite and infinite cases whilst the LS estimator's performance remains constant. This is due to the jump filtering term that absorbs most of the increments as the amplitude grows. Also, a finite-sample bias is present both estimators converge anyway as the REM performance is surprisingly consistent across paths for all estimators with 95% errors equal to a thousandth of a REM unit, even for σ = 1, 000. For a given estimator, the behaviour is similar across different topologies. The MLE performs better for sparser graphs: for σ = 1, the REM is ten times as small for the MLE on Polymer compared the LS whilst it is only twice as small on Complete.
The fact the performance of the LS estimator is not related to the underlying graph structure is difficult to interpret. We conjecture this is a side-effect of the row-normalisation applied to the dynamics matrix. However, the constant performance even for large noise amplitudes was expected as it converges even with infinite variance [25].

Robustness to changes in mesh size
The availability of data sets quantifying the same system but different mesh sizes are common (e.g. meteorological and financial data) and being able to compare the GrOU parameters across data sets with different mesh sizes becomes crucial.
Recall that the MLE converges as T N → ∞ and ∆ N → 0 as detailed in Assumption 5. The LS estimator encodes the mesh size in the estimate of exp(−Q∆ N ) as opposed to the MLE. However, for a fixed sample size N = 25, 000, we compare the estimators's performances by varying ∆ N from 0.001 to 0.200 whilst the true parameter values were generated with ∆ N = 1/12 as pictured in Figure 11.
We observe that the parameter values are linked to a particular mesh size, may it be the MLE or LS estimators as hinted by the drop in REM around the original mesh size of ∆ N = 1/12. With a larger mesh size, we see that the REM explodes exponentially (due to the log-log scale) whilst a smaller mesh size leads a larger REM which remains controlled and never exceeds 0.6.
We conclude that we ought to interpolate data to upsample the lower frequency data set rather than downsample the higher frequency one to make two sets of parameters comparable.

Conclusion
In this article, we considered the discretely-observed Graph Ornstein-Uhlenbeck process and carried out the translations of the continuous-time estimators and asymptotic results given in [21] to a non-uniform high-frequency discrete time framework. This was done by adapting the (high-frequency) double asymptotic assumptions from [36] to the multivariate context. Jump-filtered and discretised maximum likelihood estimators are shown to converge stably to the same distribution as in the continuously-observed case under a set of standard assumptions. This yields two different scales for interpretation: on a d-node graph, two-parameter θ-GrOU or d-parameter ψ-GrOU parametrisations are available; the former detailing the average behaviour of a node whilst the latter could be useful for prediction purposes. In addition, the non-uniformity of the observation grid extends the range of applicable data sets: e.g. event-based data sets closer to real-time data aggregation [47,29]. In addition to providing the consistency of the estimators, the stable convergence form a theoretical framework to perform edge pruning in conjunction to the GrOU inference [40] if the network structure is to be stochastic or time-dependent [50]. Alternatively, if the structure is deterministic but unknown, the Adaptive Lasso scheme from [21,26] extended to high-frequency observations is shown to be asymptotically normal (stably) and consistent in variable selection.
Future research comprises the assessment of the forecasting performance of such a model on high-frequency data sources (e.g. on tick-by-tick financial time series) for a large collection of systems or financial assets [29]. Note that including the impact of the neighbours' neighbours, or more generally n-th degree neighbours as in [33], could prove useful in improving the model's performance. Another direction would be the extension of the discrete-time asymptotic convergence results to Lévy-driven spatio-temporal OU processes as defined in [43] and comparing them with the coupled sparse inference on the adjacency matrix along with the GrOU inference as mentioned above. Additionally, extending the results for the studied OU-type structure for graphs to the more general continuous time autoregressive moving average (CARMA) processes [37,16] under a high-frequency sampling scheme is a promising extension. Alternatively, the application of stochastic volatility on GrOU-type process [21, Section 5] would alleviate the fixed structure of the Gaussian component and introduce empirical properties such as volatility clustering. Finally, the stable convergence results may allow to extend the GrOU process to stochastic graphs, potentially by borrowing ideas from the random graph modelling literature [50,45].

Appendix A: Lévy-Itô decomposition
In this section, we provide standard results and properties of the Lévy-Itô decomposition for Lévy processes. We say a truncation function is any R d -valued nonnegative function. The Lévy process L has the following characteristic func-tion for L at time t ≥ 0: The decomposition formulated in [13] unfolds as follows: Theorem A.1 (adapted from Theorem 3.12, [13]). Let x → τ (x) be a truncation function (e.g. τ (z) := I {x∈R d : x ≤1} (z) where z ∈ R d ). The L be a d-dimensional Lévy process with characteristic triplet (b, Σ, ν) with respect to the truncation function τ , then there exist with covariance matrix Σ and almost-surely continuous paths; such that we have uniquely for any t ≥ 0 zN t (dz) and we define the process J t := J 1 t + J 2 t . In addition, ν(dz) is a positive measure on R d \{0 d } with min( z 2 , 1)ν(dz) < ∞.
Remark A.2. In this article, we use τ (z) := I { x ≤1} (z). That is, τ is positive when z ≤ 1 allowing to capture small jumps in J 2 whilst τ is zero for larger jumps (as in the definition of b and J 1 ).
Note that N t (dz) = t 0 N (ds, dz) where N (ds, dz) is a Poisson random measure on [0, ∞) × R d . Under Assumption 6, one has L t = W t + J t . Also, the stochastic integral of J 1 is in the sense of L 0 and the second one, of J 2 t , in the L 2 sense. As noted in [20], (W t ), (J 1 t ) and (J 2 t ) are independent and the convergence in the last term is almost sure and uniform in t on compact sets. Proof of Lemma 3.1. From Lévy-Itô's decomposition, there exists a d-dimensional Brownian Motion W t with respect to P Y and covariance matrix Σ and a pure jump Lévy process J t such that L t = W t +J t (without drift by assumption). The continuous P 0 -martingale part is then expressed as Recall that ψ = vec(Q) and that Since this is true for any i, j ∈ {1, . . . , d}, we conclude that Using Equation (B.1.1), we obtain: where s N := t k such that t k ≤ s < t k+1 for k ∈ {1, . . . , N }. Therefore, we write respectively the terms above as As in the univariate case [36, Section 3, p. 929], we will prove that Z 1 , by stationarity.
We prove that K N approximates its continuously-observed counterpart (K T N ) in the L 1 sense as follows: , by stationarity, Again, using both the triangle and Cauchy-Schwarz's inequalities and finally the Lévy-Itô decomposition, we see that we have s. as t → ∞ by ergodicity and therefore Next, by the continuous mapping theorem, It remains to show that Z 2 which is proved in a similar fashion albeit using a multivariate martingale central limit theorem, namely Theorem 2.2, [22]. Again, we have that We shall now focus on the remainder term T Similarly, we find the L 1 convergence of those quantities as N → ∞. For i, j, k, l ∈ {1, . . . , d}, we have that (similarly as above) This is justified by the ergodicity of (Y t , t ≥ 0) and path continuity of (M t , t ≥ 0). Therefore, by Theorem 2.2, [22], U t ·M t converges F-stably to N 0 We obtain the F-stable convergence by Slutsky's lemma: which concludes the proof.
Appendix C: Proof for the Adaptive Lasso Proof. We first prove the asymptotic normality of the estimator before proving the property of consistency in variable selection. The likelihood part of the objective function is to be written as We then reformulate further the likelihood. In the proof of Lemma 3.1 (Appendix B), we show that . We recall that Q N − Q N P − − → 0 d 2 and that, from that same we also obtain Using those notations, the log-likelihood is reformulated as We then proceed similarly to the proof of Th. 5.3.1, [21] (see App. A.10 therein). First, without loss of generality, we change the penalty rate from λ to λT N since it does not depend on Q. Then, by writing where the third term κ 3 N represents the discretisation error which will become asymptotically negligible as shown next. We denote and define the complement of Q 0 as follows Q 0 := {(k, l) ∈ Q 0 : 1 ≤ k, l ≤ d}. Then, we treat each term κ i N , i ∈ {1, 2, 3}, separately: • For the first function κ 1 N : by Slutsky's lemma (for stable convergence [28, Th. 1.1]), the first term converges stably to a standard (centred) Gaussian distribution with covariance K 1 [21] and [26] use a d × d Gaussian matrix Z with mean zero and covariance Cov(vec(Z), vec kl . Therefore, for any M ∈ M d (R), tr(MZ) is a Gaussian random variable with mean zero and variance Also, the second term converges almost-surely to −K 1 ∞ /2 such that by Slutsky's lemma, we have for any • For the second function κ 2 N : the argument from Appendix A.10, [21] can be replicated exactly and we do not present it for brevity. In essence, by consistency of the MLE (Corollary 3.2) and since λ(N )N 1/2 → 0 and λ(N )N (γ+1)/2 → ∞ as N → ∞, the objective function diverges to −∞ if supp(M) ∩ G 0 = ∅ and zero otherwise for any M ∈ M d (R).
• For the third function κ 3 N : each term converges to zero in probability for any M ∈ M d (R) as N → ∞.
We summarise the resulting asymptotic properties as follows: This is the same asymptotic property is in [21] and therefore, we borrow their conclusion to show that the maximum in the limit t → ∞, denoted M, is given by Given the covariance structure of Z, vec( M |G0 ) is a centred Gaussian random vector with covariance (E(Y ∞ Y ∞ ) −1 ⊗ Σ) |G0×G0 which concludes the proof of asymptotic normality. The consistency in variable selection is shown similarly to the continuoustime observations case [21]: the asymptotic normality of Q AL,N on G 0 yields that P ( Q AL,N ) ij = 0 → 1 if (i, j) ∈ G 0 . We then show that P ( Q AL,N ) ij = 0 → 1 if (i, j) ∈ G 0 . We suppose that ( Q AL,N ) ij = 0 and then, by taking the derivative of the objective function with respect to the (i, j)-th parameter of Q AL,N .
We multiply by T −1/2 N , set it to zero before taking the absolute value on both sides. We obtain that For the right-hand side: from the first part of the proof, we have T For the left-hand side: since s., the first term is normally-distributed with finite variance as N → ∞. By the asymptotic normality of T 1/2 N ( Q AL,N ) |G0 , the second term is also asymptotically normal with finite variance. That is, the left-hand side is the absolute value of the sum of two Gaussian random variables whose probability to be larger or equal to the right-hand side (which diverges to ∞) tends to zero with probability one. The right-hand side was computed under the assumption that ( Q AL,t ) ij = 0, and the probability that event happening is upper-bounded by zero, hence we have that P ( Q AL,t ) ij = 0 → 1 if (i, j) ∈ G 0 . This concludes the proof of the consistency in variable selection.

Appendix D: Technical results
In this section, we present a collection of intermediate results leading up to Theorems 3.3-3.4 & 3.8-3.9 whilst the proofs are presented in Appendix E and F, respectively, for the finite and infinite jump activity case.

D.1. Intuition
From Lemma 3.1 in Section 2.4.2, the discretised unfiltered estimator ψ N converges (stably) in distribution to the limit of the continuous-time estimator with time horizon T N , namely ψ T N , in the sense that The purpose of the following sections is to prove that we can extend this result to the jump-filtered estimator ψ N . We show that filtering the jumps with a particular threshold gives way to a consistent estimation in both the finite and infinite jump activity cases. This is carried out by treating the different terms given by the Lévy-Itô decomposition of L separately. Note that the decomposition differs significantly in terms of the components' properties be it in the finite or the infinite jump activity case. The latter requires intermediate filtering steps to relate the discretised jump-filtered estimator and its unfiltered counterpart.

D.2. Finite jump activity
In this section, we prove Theorems 3.3 and 3.4. Recall that, under Assumption 6, the Lévy-Itô decomposition (see Appendix A) yields that there exists a centred Gaussian Lévy process W t with covariance matrix Σ and almost-surely continuous paths as well as a pure-jump Lévy process J t independent of W t such that Finally, under the finite jump activity assumption (Assumption 8), remark that J t is a compound Poisson process. Indeed, there are Poisson processes N For simplicity, we introduce some notation for the drift term: Notation D.2. We denote by ∆ k D the drift term between t k and t k+1 as expressed by We first introduce results proved in [36] but adapted to multiple dimensions. A crucial result is the upper bound on the probability of the continuous part increments {|∆ k W (j) + ∆ k D (j) | > ∆ 1/2−δ N } for any δ ∈ (0, 1/2) as formulated below: Lemma D.3. (adapted from Lemma 3.7, [36]) Suppose that for any i ∈ {1, . . . , d}, Then, for any δ ∈ (0, 1/2) and k ∈ {1, . . . , N }, we have: Proof. The proof in presented on page 925, [36].
[36] introduced a family of events when a consensus between the filtering and the absence of jumps takes place. We extend to our multivariate framework by introducing equivalent events componentwise.
Definition D.4. Let (Y, t ≥ 0) be a GrOU process such that Assumptions 1-8 hold. We define the sequence of events (A (i) k,N : k ∈ {1, . . . , N − 1}) for the i-th node where a small increments corresponds to the absence of jumps, i.e. define Notation D.5. We sometimes write I{S} for ω → I S (ω) where S ∈ F.
They proved that the filtering and the absence of jumps coincide all the time with probability one as N → ∞ in the following sense: Proof. See Section 3.3.1, pages 925-927, [36].
We now extend Lemma 3.9, [36] to the multivariate setting. Let i, j ∈ {1, . . . , d}. We show that we can approximate (in the L 1 sense) the discretised integral of Y (i) with respect to the continuous part of Y N , then we have as N → ∞: for any i, j ∈ {1, . . . , d}.
Proof of Lemma D.7. Let i, j ∈ {1, . . . , d}. Following Lemma D.6, by conditioning on A (j) k,N , i.e. where increments smaller than the jump thresholds coincide with the absence of jumps at all times, we have: happens. Thus, we define the event B (j) k,N := {∆ k N (j) > 0}. Now, we can rewrite the above equality as follows By the Lévy-Itô decomposition, the increments of the continuous part can be written Hence, by the triangle inequality Thus, we compute both terms of the right-hand side separately. As in the proof of Lemma 3.9, [36], observe that since ∆ k N (j) is Poisson-distributed with rate as well as between ∆ k N (j) and Y (i) k gives way to a succinct decomposition For the second term, we use Hölder's inequality to obtain Then, by independence of Y (i) k and ∆ k N (j) , we have: Again by Hölder's inequality, observe that: Since Y is stationary and has finite second moments, we obtain by Fubini's theorem We conclude that Lemma D.7 is essential as it links the continuous part of Y increments to the jump-filtered increments. We can then compute the difference ψ N − ψ N element-wise and show that as given in Appendix E.1:

D.3. Infinite jump activity
Similarly to Section 4, [36], we put forth a collection of convergence results to prove Theorems 3.8 & 3.9: we first justify the negligibility of some jump-filtered quantities (in the sense of convergence in probability to zero). We also prove that one can properly estimate various unfiltered quantities with jump-filtered equivalents with a high-frequency sampling scheme. Since we are to show the consistency of the jump-filtered estimator with respect to the unfiltered one by splitting the difference in four different terms, the proofs are relegated to Appendix F for clarity.
The jump-filtered part of the process is actually going to zero in probability as shown in the following lemma: Lemma D.8. (Extension of Lemma 4.14, [36]) Under the assumptions of Theorem 3.8, we have for i, j ∈ {1, . . . , d}: Proof. See Appendix F.1.
The second lemma proves that filtered signals with large jumps are negligible: Lemma D.9. (Extension of Lemma 4.13, [36]) Under the assumptions of Theorem 3.8, we have for i, j ∈ {1, . . . , d}: Finally, the last two results of this section allow to asymptotically filter out the contribution of small jumps larger than v (j) n /2, for j ∈ {1, . . . , d}. First, one recalls an essential couple of results from [36]: Lemma D.10. (Equation (38), [36]) Under the assumptions of Theorem 3.8 and for any j ∈ {1, . . . , d}, one obtains Also, one has that P v (j) Proof. See page 949, [36].
Secondly, we prove the negligibility of a J 2 −filtered discrete equivalent to A N (Definition 2.3) as a convergence in probability to zero in the following sense: Lemma D.11. (Extension of Lemma 4.15. [36]) Under the assumption of Theorem 3.8 and for i, j ∈ {1, . . . , d}, we have: It can be shown that in infinite jump activity case, small increments of Y are dominated by small jumps. We prove that sums of jumps increments are asymptotically zero if one filters for small enough jumps: Lemma D.12. Under the assumptions of Theorem 3.8, we have for i, j ∈ {1, . . . , d}: Proof. See Appendix F.4.
As an intermediary result, we prove that one can focus on the small-jump component as non-zero large jumps are asymptotically negligible as follows Lemma D. 13. Under the assumptions of Theorem 3.8, we have for i, j ∈ {1, . . . , d}: Proof. See Appendix F.5.
Then, we prove that filtering for small increments is equivalent to filtering for small jumps in the sense of the following lemma: Lemma D.14. (Extension of Lemma 4.12, [36]) Under the assumptions of Theorem 3.8, we have for i, j ∈ {1, . . . , d}: Also, we prove that more generally Proof. See Appendix F.6.

D.3.1. Filtering the drift
The following lemma shows that under jump-filtering with threshold v (j) N , the discrete integration with respect to the drift can be approximated regardless if the filtering is applied.

D.3.2. Filtering the continuous martingale component
To show that filtering out jumps is a valid option in the framework of Theorem 3.8, we propose the lemma that follows: Lemma D.16. (Adapted from Lemma 4.10, [36]) In the framework of Theorem 3.8, for i, j ∈ {1, . . . , d} and as N → ∞: Proof. See Appendix F.8.
We have provided two collections of lemmas laying out the proofs for the central limit theorems for both the finite and infinite jump activities. They consist of incremental filtering and decompositions gradually linking the unfiltered to the filtered estimators. In the following section, we prove that the discretised unfiltered estimator is consistent as presented in Lemma 3.1.

Appendix E: Proofs for the finite activity case
In this technical section, we provide the proofs of the finite activity intermediate results presented in Appendix D.2.

E.1. Proof of Lemma 3.6
Proof of Lemma 3.6. Observe that According to Lemma D.7, for any i, j ∈ {1, . . . , d}, we have under P 0 Also, by ergodicity and the continuous mapping theorem, one has that s. which is finite componentwise and this concludes the proof.

Appendix F: Proofs for the infinite activity case
In this technical section, we provide the proofs of the infinite activity intermediate results presented in Appendix D.3.

F.1. Proof of Lemma D.8
To prove this result, we introduce some new notation and define for X ∈ {Y, W, D, J 1 , J 2 } and i, j ∈ {1, . . . , d}. We prove the convergence in probability to zero of B X N (i, j) for each X ∈ {Y, W, D, J 1 , J 2 }, starting with B W N (i, j) and B D N (i, j) in the section below.

F.1.1. Intermediate results
We prove two results on the limiting behaviour of, respectively, B W N and B D N as follows Lemma F.1. Under the assumptions of Lemma D.8, we have for i, j ∈ {1, . . . , d} Proof of Lemma F.1. B W N can be bounded using Hölder's inequality and secondly by independence between Y (j) k and ∆ k W (j) : , by Lemma D.3 with δ = 1/2 − β (j) and l (j) = 2, As we bound the L 1 -norm of B W N (i, j), we obtain the convergence in probability to zero by Assumption 7.
Lemma F.2. Under the assumptions of Lemma D.8, we have for i, j ∈ {1, . . . , d} Proof of F.2. We prove this result similarly to Lemma F.1. First, we recreate increments as follows: The first double sum is upper-bounded by Hölder's inequality and Fubini's theorem as follows: Hence, we have: The second term can be processed in a similar fashion. Using Hölder's inequality and Assumption 9 (i), we obtain: Therefore, we have that: Hence, as we bound the L 1 -norm of B D N (i, j), we obtain the convergence in probability to zero by Assumption 7.
N . By the triangle inequality, we then split the quantity of interest First, Lemmas F.1 and F.2 yield that Finally, remark that by independence and again by Hölder's inequality, B J 1 N (i, j) and B J 2 N (i, j) can be treated similarly to B W N (i, j). All four part are then dominated by T

F.2. Proof of Lemma D.9
Proof of Lemma D.9. The key idea of the proof is to prove that the indicator by adding and subtracting ∆ k Y (j) in the first term. That implies that, as given in [36]: Indeed, if that were not the case, then one would have N leading to a contradiction in Equation (F.2.1) when taking the absolute value. Therefore, and Lemma D.8 allows to conclude directly.
For the first term on the right-hand side, we repeat a similar argument to the proof of Lemma F.2, presented in F.1: by Assumption 9-(i) and Hölder's inequality, we obtain that Next, for the second term, the independence between Y (i) k and ∆ k J 2,(j) yields that ). The fourth element U J 2 N (i, j) is bounded using Assumption 9-(iii) in the sense that for some 0 > 0, for N large enough v (j) Then, the independence between Y Also, recall that Lemma D. 10 gives ). Therefore, we obtain: Hence, by Markov's inequality and since α ∈ (0, 2), U J 2 N (i, j) converges to zero in probability as N → ∞.
Finally, for the third term U J 1 N we use the independence between Y (i) k and ∆ k J 1,(j) . By definition E |∆ k J 1,(i) | = O(∆ N ), thus: as N → ∞, all four terms converge to zero in L 1 hence in probability as N → ∞ which concludes the proof.
Using Hölder's inequality and the independence between Y (i) k , ∆ k J 1,(j) and ∆ k J 2,(j) , the part with respect to ∆ k J 1,(j) is upper-bounded as follows: where we used the assumption N ∆ N = O(T N ) as N → ∞. For the term with respect to ∆ k J 2 , we proceed using Lemma D.10: which concludes the proof.

F.5. Proof of Lemma D.13
Proof of Lemma D. 13. We prove that the pure jump part J = J 1 + J 2 is dominated by small jumps in the sense that the J 1 -jump contribution is asymptotically zero for filtered jumps as follows Again, using ∆ k Y (j) = ∆ k D (j) + ∆ k W (j) + ∆ k J (j) , one has: As remarked on page 943, [36], J 1,(j) is a compound Poisson with a finite intensity given by ν(R\[−1, 1]) < ∞ such that one can write P |∆ k J 1,(j) | = 0 = O(∆ N ). This implies that The second item can be processed in a similar way to the term U D N (i, j) in the proof of Lemma D.8 (as in Lemma F.2, Appendix F.1). Hence since P |∆ k J 1,(j) | = 0 = O(∆ N ). Finally, by Lemma D.12, the third term is asymptotically zero w.r.t. to the L 1 -norm as N → ∞. Therefore: F.6. Proof of Lemma D.14 Proof of Lemma D.14. Consider the whole pure jump process J (j) instead of J 2,(j) . We split this quantity of interest into two parts as follows: We note that the first term is exactly the quantity converging in probability to zero from Lemma D.9 as N → ∞. By Lemma D.13, we know that (F.6.1) Next, we show that Indeed, observe that if ∆J 1,(j) = 0, then ∆ k J (j) = ∆ k J 2,(j) ; thus Similarly to the proof of Lemma 4.12, [36], we split the event in the indicator function as follows: which redirects the convergence proof to the results of Lemmas D.8 and D.11. This inclusion is true since either we have |∆ k J 2,(j) | < v (j) N /2 and in this case N . Given this convergence result and (F.6.1), we obtain that which proves the desired result.

F.7. Proof of Lemma D.15
Proof of Lemma D.15. Remark that Lemma D.14 yields that we can prove the result for the quantity above where Similarly to the proof of Lemma D.8 (Lemma F.2, presented in Appendix F.1), we write the following up to a term that converges to zero in probability as N → ∞: Since v (j) ) and by Chebychev's inequality, we have P |∆ k J (j) | > 2v and ∆ k J (j) , we obtain directly for the second term: and then for the first term, according to Hölder's inequality, we obtain  Proof of Lemma D. 16. We define the process containing all components Y except the infinite-activity pure jump process J 2 , namely define Y t := Y 0 − t 0 QY s ds + W t + J 1 t . Consider splitting Y t into Y t − Y t and J 2 t as follows: The first term V 1 N (i, j) is similar quantity to Lemma D.7 and by definition V 1 N (i, j) has finite activity jump. By applying this result, we have that V 1 N (i, j) = O p (∆ 1/2 N T N ): it converges in probability to zero as N → ∞. The second term can be rewritten if one remarks that I{|∆ k Y (j) | ≤ v (j) with obvious definitions for {V 2,l N (i, j), : l ∈ {1, 2}}. The first term V 2,1 N (i, j) can be shown to converge to zero in probability since, informally, it implies that the process J 2 is negligible overall as below v (j) N = ∆ β (j) N at every time step. Indeed, by sub-additivity: Since we now have a univariate right-hand side, we now use the proof of Lemma 4.10, [36]. Remark that by the triangle inequality Recall that, by Lemma D.3 with δ = 1/2 − β (j) and l (j) = 2, we have ).
Since this sum goes to zero as N → ∞, when a large increment in Y occurs, it is highly likely that a large jump occurred (i.e. |∆ k J 1,(j) | = 0). Also, again by the N }. Using the independence between ∆ k J 1,(j) and ∆ k J 2,(j) , we prove that V 2,1 N (i, j) is negligible as follows: ) is given by the Chebychev's inequality. The term V 2,2 N (i, j) can undergo the same treatment and we have Again, we fall back to the univariate case and prove the convergence as presented in [36]. Recall that J 1 is a compound Poisson process with (N (i) t , i ∈ {1, . . . , d}) as a collection of d Poisson processes with counting processes (N (i) : i ∈ {1, . . . , d}). According to (F.8.1), we have that Since J 1 is defined as the large-jump component, we show that it has a negligible impact when {|∆ k Y (j) | ≤ 2v which converges to zero in the limit. This shows that we have ∆ k Y (j) = ∆ k W (j) + ∆ k D (j) a.s. on {|∆ k Y (j) | ≤ 2v (j) Next, according to Lemma D.15 (and its proof), the part of V 2,2 N (i, j) with respect to ∆ k D (j) is such that The second part of V 2,2 N (i, j) is given by Recall the independence between Y (i) k , ∆ k W (j) and ∆ k J 2,(j) . Its expectation is then zero and its L 2 -norm we can bound as follows by : Therefore, we have also proved that this second term of V 2,2 N (i, j) converges to zero in probability and, as a result, so does V 2,2 N (i, j) and, in turn, V 2 N (i, j) as N → ∞.
Finally, for the third term, recall Lemma D.13 such that we consider Then, by the proof of Lemma D.14, we know that Since v (j) N → 0 as N → ∞, then by Assumption 9-(iii), we have that E ∆ k J 2,(j) I{|∆ k J 2,(j) | ≤ 2v (j) N } = 0, for N large enough.
Then, by independence of Y (i) k and ∆ k J 2,(j) , we obtain that Finally, again by independence and by the triangle inequality, the second moment of this quantity is bounded by This concludes the proof.
F.9. Proof of Lemma 3.7 Proof of Lemma 3.7. Observe that Now, by Lemma D. 16, We know that for any i, j ∈ {1, . . . , d}, we have that Also, recall that S N := K N ⊗I d×d where K N = which is finite componentwise. Therefore, similarly to Lemma 3.6, we obtain that T 1/2 F.10. Proof of Theorem 3.8 We can finally prove Theorem 3.8 using the lemmas above.
Proof of Theorem 3.8. By Lemma 3.1, we know that the discretised unfiltered estimator ψ N converges in distribution as follows: By Lemma 3.7 we have that We conclude by Slutsky's lemma.