Valentin Blot
email: valentin.blot@inria.fr

Gilles Dowek
email: gilles.dowek@ens-paris-saclay.fr

Thomas Traversié
email: thomas.traversie@inria.fr

An Implementation of Set Theory with Pointed Graphs in Dedukti

DEDUKTI is a type-checker for the λ Π-calculus modulo theory, a logical framework that allows the extension of conversion with user-defined rewrite rules. In this paper, we present the implementation of a version of Dowek-Miquel's intuitionistic set theory in DEDUKTI. To do so, we adapt this theory -based on the concept of pointed graphs -from Deduction modulo theory to λ Π-calculus modulo theory, and we formally write the proofs in DEDUKTI. In particular, this implementation requires the definition of a deep embedding of a certain class of formulas, as well as its interpretation in the theory.

Introduction

During the last decades, theorem provers have attracted a lot of interest from various fields of science. Especially, the use of such tools has become more and more common in the areas of software verification and formalization of mathematics. As an example, verification has become mandatory for software running on aircrafts, and formalization of mathematical results allowed the scientific community to trust some complicated proofs of recent results, as well as identify and fix errors in some cases.

This growing interest has triggered the development of many theorem provers, with various focuses and based on a large range of theories. While the diversity of theorem provers provides users with a wide range of tools, this comes at the expense of a lack of interaction and reusability between proof fragments done in different tools. The λ Π-calculus modulo theory is a logical framework supporting dependent types and user-definable rewrite rules, and DEDUKTI [START_REF] Assaf | Dedukti: a Logical Framework based on the λ Π-Calculus Modulo Theory[END_REF] is a type-checker for this framework. The definition of well-chosen rewrite rules allows the encoding of various logics, and hence the translations between various theorem provers through the DEDUKTI tool.

Several theorem provers, such as MIZAR, ATELIER B and ISABELLE/ZF, are based on set theory. In order to extend the interoperability allowed by DEDUKTI to these provers, it is necessary to encode set theory in the λ Π-calculus modulo theory, and to implement this encoding in DEDUKTI. The goal of this paper is to present such an encoding and implementation. In order to facilitate the proof of its correctness, we implement this encoding in the tool LAMBDAPI that provides tactics to help the user in the production of proofs type-checkable in DEDUKTI.

Stating each axiom of set theory in DEDUKTI would lead to an implementation that does not satisfy a cut-elimination theorem. In particular, it would forbid extraction of witnesses from constructive existence proofs.

An alternative option would be to orient these axioms as rewriting rules. For instance the powerset axiom x ∈ P(y) ⇔ x ⊆ y would be replaced with the rewrite rule x ∈ P(y) -→ x ⊆ y. However, as pointed out by Crabbé [START_REF] Crabbé | Non-normalisation de la théorie de Zermelo[END_REF], such a formulation of set theory does not satisfy a cut-elimination property either. Indeed, if for a given set a we define the set b ≡ {x ∈ a | x ∈ x} then b ∈ b would get rewritten to b ∈ a ∧ b ∈ b, leading to an infinite reduction.

In the current paper, we represent sets as pointed graphs [START_REF] Aczel | Non well-founded sets[END_REF]. With such a formulation, we can prove a proof normalization theorem: for every proof in natural deduction there exists a normal proof of the same statement. Such an encoding of set theory has been defined in the context of Deduction modulo theory [START_REF] Dowek | Cut elimination for Zermelo set theory[END_REF], together with pencil and paper proofs. At that time, the λ Π-calculus modulo theory had not been defined yet, and the first contribution of the current paper is to adapt this encoding to the λ Π-calculus modulo theory. In particular, we avoid the original definition of classes of nodes by using the quantification on propositions permitted by λ Π-calculus modulo theory. This reduces the size of the signature from 31 to 26 symbols.

The second contribution consists of the implementation of this new theory in DEDUKTI. In the original formulation of the theory some lemmas are only valid for a specific class of formulas, as their proofs proceed by induction on the structure of these formulas. Implementing these lemmas in DEDUKTI requires the definition of an inductive sort of formulas together with an interpretation of these formulas into the general type of propositions. This interpretation is defined with rewriting rules. The current development is the first, to our knowledge, to use such a reflection principle defined with rewriting rules in DEDUKTI. The generality of this method still remains to be investigated.

The Theory of Pointed Graphs

The set theory with pointed graphs developed by Dowek and Miquel is called Intuitionistic Zermelo set theory in Deduction modulo (IZmod). We give here an informal presentation of the ideas developed in [START_REF] Dowek | Cut elimination for Zermelo set theory[END_REF].

Sets as Pointed Graphs

In the IZmod theory, sets are represented by pointed graphs, that is, directed graphs with a distinguished node: the root. We give here the representation of several sets as pointed graphs. The root is indicated by the filled circle. The location of the root is important as it changes the set represented by the pointed graph. For example, in the third graph, the node b is irrelevant since the root is c. Distinct pointed graphs can represent the same set. For example, here are various representations of ordinals 0 = / 0, 1 = { / 0} and 2 = { / 0, { / 0}} as pointed graphs. IZmod has a sort of pointed graphs and a sort of nodes. If a is a pointed graph, root(a) is the root of a, which is a node. If x is a node of a, then a/x refers to the pointed graph a where x is the new root (this is a re-rooting operation). If y is also a node of a, x η a y is the proposition asserting that there is a edge in a from y to x. Following this interpretation, IZmod contains the following rewriting rules:

x η a/z y -→ x η a y root(a/x) -→ x (a/x)/y -→ a/y

As noted above, distinct pointed graphs can represent the same set. Two pointed graphs representing the same set are related through a notion of bisimilarity, denoted ≃, and defined with the rewrite rule:

a ≃ b -→ ∃r, r root(a) root(b) ∧ ∀x∀x ′ ∀y (x ′ η a x ∧ r x y ⇒ ∃y ′ (y ′ η b y ∧ r x ′ y ′ )) ∧ ∀y∀y ′ ∀x (y ′ η b y ∧ r x y ⇒ ∃x ′ (x ′ η a x ∧ r x ′ y ′ ))
As an example, the following pointed graphs are bisimilar, hence they represent the same set: In set theory, there only exists one sort: the sets. The IZmod theory has two main sorts -nodes and pointed graphs -and two notions related to set membership -relations between pointed graphs and relations between nodes within the same graph. As an example, take / 0 ∈ { / 0}:

k ℓ ∈ k ℓ
Since there is an edge from k to ℓ, the set represented by the pointed graph with root ℓ is an element of the set represented by the pointed graph with root k. But any pointed graph bisimilar to the pointed graph with root ℓ also represents a set that is an element of the set represented by the pointed graph with root k. This leads to the definition of a membership relation ∈ with the rewrite rule:

a ∈ b -→ ∃x (x η b root(b) ∧ a ≃ (b/x))
In order to define set constructions such as pairs or powersets, we need a way of joining two existing graphs together, as well as a way of adding a new node to an existing graph. To this end, the IZmod theory has two disjoint injections i and j on nodes, and a symbol for some constant node o that is not in the image of i or j. There are also inverses i ′ and j ′ to i and j, and predicates I and J on the images of i and j:

i ′ (i(x)) -→ x I(i(x)) -→ ⊤ I( j(x)) -→ ⊥ I(o) -→ ⊥ j ′ ( j(x)) -→ x J( j(x)) -→ ⊤ J(i(x)) -→ ⊥ J(o) -→ ⊥
For example, we want to define a pointed graph {a, b} that satisfies the axiom of pairing, that is to have ∀x (x ∈ {a, b} ⇔ (x ≃ a ∨ x ≃ b)). The intuitive idea to build this pointed graph is to create a new node that will be parent of the root of a and of the root of b, and to define this new node to be the root of {a, b}. We need to ensure that all the nodes of {a, b} are different from each other. To comply with this constraint, we choose o for the new node and we use injection i on the nodes of a and injection j on the nodes of b. This construction is formalized by the rewriting rules

root({a, b}) -→ o x η {a,b} x ′ -→ (∃y∃y ′ (x = i(y) ∧ x ′ = i(y ′ ) ∧ y η a y ′ )) ∨ (∃y∃y ′ (x = j(y) ∧ x ′ = j(y ′ ) ∧ y η b y ′ )) ∨ (x = i(root(a)) ∧ x ′ = o) ∨ (x = j(root(b)) ∧ x ′ = o)
We can proceed similarly for other constructions such that powerset or union.

Set Theories

The membership ∈ and bisimilarity ≃ relations define a set theory IZst that lies between Zermelo (Z) and Zermelo-Fraenkel (ZF). This theory does not include the Replacement scheme but it contains two additional axioms: Strong Extensionality and Transitive Closure:

Strong Extensionality axiom. ∀x 1 ...∀x n ∀a∀b (R(a, b) ∧ ∀x∀x ′ ∀y (x ′ ∈ x ∧ R(x, y) ⇒ ∃y ′ (y ′ ∈ y ∧ R(x ′ , y ′ ))) ∧ ∀y∀y ′ ∀x (y ′ ∈ y ∧ R(x, y) ⇒ ∃x ′ (x ′ ∈ x ∧ R(x ′ , y ′ ))) ⇒ a ≃ b)
where R(a, b) is a formula with free variables x 1 , ..., x n . Note that the hypothesis of the Strong Extensionality axiom mimics the structure of the rewrite rule for ≃.

Transitive Closure axiom. ∀a∃e (a ⊆ e ∧ ∀x∀y (x ∈ y ∧ y ∈ e ⇒ x ∈ e))
The Transitive Closure axiom conveys the idea that every set is included in a transitive set.

Strong Extensionality can be deduced from Foundation and Transitive Closure can be derived from Replacement. Moreover, Strong Extensionality implies Extensionality 1 . Therefore we have the following picture:

Z

Pairing Union Extensionality

Powerset Comprehension Infinity

IZst

Strong Extensionality Transitive Closure

ZF

Foundation Replacement

The IZmod theory is an extension of IZst set theory with an encoding of pointed graphs. Indeed, the axioms of IZst are provable in IZmod [5, see Section 3.3].

The Language of Pointed Graphs 3.1 Sorts

There are four sorts in IZmod [5, see Section 3.2]: graphs, nodes, classes of nodes and binary relations between nodes. We implemented the first two via the standard technique [START_REF] Blanqui | Some Axioms for Mathematics[END_REF] of defining a universe Set : TYPE of sorts, a function El of type Set → TYPE and two constants graph and node of type Set: constant symbol graph : Set ; constant symbol node : Set ;

1 See Section 5. [START_REF] Blanqui | Some Axioms for Mathematics[END_REF] We also introduce a constant arrow of type Set → Set → Set, together with its associated rewrite rule El (x arrow y) -→ (El x) → (El y).

Contrary to Deduction modulo, classes of nodes and binary relations on nodes can be expressed primitively via quantifications on propositions in λ Π-calculus modulo theory. The sort of classes is defined in DEDUKTI as El node → El prop and that of binary relations as El node → El node → El prop.

The symbols graph and node are specific to the expression of IZmod in DEDUKTI, while in contrast the symbols Set, El and prop are part of the standard library of DEDUKTI.

Signature

The signature of IZmod contains 33 symbols (the 31 from [5, see Table 2], plus the symbol of empty set [5, see Section 3.2] and the symbol of inductive set [5, see Section 2.1]). As we have replaced the sorts for classes and relations with primitive DEDUKTI types, we do not need specific predicate symbols mem and rel to apply a class to a node or a relations to two. By virtue of the encoding of classes of nodes as elements of El node → El prop, mem(x, P) can simply be expressed as P x, and similarly rel(x, y, r) can simply be expressed as r x y. In the same way, we do not need symbols g x,y 1 ,...,y n ,P and g ′

x,x ′ ,y 1 ,...,y n ,P to build classes and relations. Finally, the equality symbol is part of the standard library of DEDUKTI. The signature is thus reduced to 28 symbols. Below are the first 25 symbols, the specific case of the comprehension symbol is treated later. Some symbols are defined as constant since they are not subject to rewriting rules.

Rewriting Rules

Among the rewrite rules [5, see Table 3], we can drop the compatibility of equality with mem as it comes for free with the Leibniz equality of DEDUKTI's standard library. We can also drop the rewrite rules determining the behavior of g x,y 1 ,...,y n ,P and g ′ x,x ′ ,y 1 ,...,y n ,P , as they are superseeded by our encoding of classes and binary relations. The remaining rewrite rules are easy to implement in DEDUKTI: General rule eta ( cr $a $z ) $x $y → eta $a $x $y ; rule root ( cr $a $x ) → $x ; rule ( cr ( cr $a $x ) $y ) → cr $a $y ; The pointed graph omega represents the set of von Neumann ordinals [5, see Table 2]. Indeed, the rewriting rule of eta omega x x' states that in omega the edges correspond to the order relation < and that there is an edge from every natural number to the root of omega. Therefore, the elements of omega are the natural numbers, and the elements of some natural number are all the natural numbers strictly below it, which corresponds to the definition of von Neumann ordinals.

Relocations rule i ' ( i $x ) → $x ; rule j ' ( j $x ) → $x ; rule ρ ' (ρ $x ) → $x ; rule I ( i $x ) → ⊤; rule J ( j $x ) → ⊤; rule I ( j $x ) → ⊥; rule J ( i $x ) → ⊥; rule I ( o ) → ⊥; rule J ( o ) → ⊥; rule Pred ( Succ $x ) → $x ; rule Null O → ⊤; rule Nat O → ⊤; rule Null ( Succ $x ) → ⊥; rule Nat ( Succ $x ) → Nat $x ; rule $x < O → ⊥; rule $x < ( Succ $y ) → ( $x < $y ) ∨ ( $x = $y );

The Language of Formulas

This section is devoted to the operator that builds sets by comprehension. Indeed, in IZmod, the only valid sets of the form {x ∈ a | P(x)} are those for which P does not contain elements specific to the implementation of sets as pointed graphs, that is, all quantifiers in P must be on graph and P can only use symbols ∈, ≃ and the logical connectives [5, see Table 5]. Taking this restriction into account is the main challenge of our implementation. In order to achieve this goal, we define a datatype representing the class of formulas just described, as well as an interpretation of this datatype into El prop.

Formulas

We define the constant f ormula of type Set and the operators for building elements of El f ormula. Variables are represented as natural numbers and we do not need to handle α-conversion as we do not need substitution. We also define an induction principle over first-order formulas: constant symbol recF : Π ( P : El formula → El prop ) , π ( '∀ x , '∀ y , P ( eqF x y )) → π ( '∀ x , '∀ y , P ( inF x y )) → π ( '∀ f , '∀ g , ( P f ∧ P g )

⇒ ( P ( andF f g ))) → π ( '∀ f , '∀ g , ( P f ∧ P g ) ⇒ ( P ( orF f g ))) → π ( '∀ f , '∀ g , ( P f ∧ P g ) ⇒ ( P ( impF f g ))) → π ( '∀ f , ( P f ) ⇒ ( '∀ x , P ( allF x f ))) → π ( '∀ f , ( P f ) ⇒ ( '∀ x , P ( exF x f ))) → π ( P tF ) → π ( P fF ) → π ( '∀ f , P f );
The symbol π is used to map each proposition (of sort El prop) to the type of its proofs (of sort TYPE).

We can finally prove that interpretation is invariant under ≃-equivalent valuations:

Theorem 4.2 ∀ f ∀σ ∀σ ′ ((interpretation σ f ∧ ∀x (σ x ≃ σ ′ x)) ⇒ interpretation σ ′ f )
Proof. By induction over formula f . □

Comprehension, Empty Set and Inductive Set

With this syntax for restricted formulas at hand together with the interpretation in prop, we are finally able to define the comprehension construction in DEDUKTI:

symbol comp : El graph → ( El nat → El graph ) → El formula → El graph ;
together with its rewrite rules

rule eta ( comp $a $σ $f ) $x $x ' → ( '∃ y , '∃ y ' , (( $x = i y ) ∧ ( $x ' = i y ') ∧ eta $a y y ')) ∨ ( '∃ y , ( $x = i y ) ∧ ( $x ' = o ) ∧ ( eta $a y ( root $a )) ∧ ( interpretation ( update $σ zero ( cr $a y )) $f )); rule root ( comp $a $σ $f ) → o ;
Moreover, we are now able to define two symbols related to the lemmas of the Infinity section [5, see Table 5]: the empty set symbol and the symbol for inductive sets Ind [5, see Section 2.1]. We implement empty_set of type El graph and Ind of type El graph → El prop. We use comprehension with the set of natural numbers and the formula f F to define the rewrite rule of the empty set: rule empty_set → comp omega (λ _ , omega ) fF ;

We can now define inductive sets:

rule Ind $c → ( empty_set ∈ $c )
∧ ( '∀ a , ( a ∈ $c ) ⇒ (( join ( pair a ( pair a a ))) ∈ $c ));

Inductive sets are used in the axiom of infinity. The previous rule states that c is inductive if (i) it contains the empty set and (ii) if a ∈ c then a ∪ {a} ∈ c. Here we use the fact that ∪{a, {a, a}} = ∪{a, {a}} = a ∪ {a}.

Lemmas

Dowek and Miquel proved that the theory IZmod of pointed graphs does validate all the axioms of IZst. This required 53 lemmas that were informally proved: each axiom of IZst corresponds to a lemma in IZmod, and we have additional intermediary lemmas on the structure of pointed graphs [5, see Tables 4 and5]. For instance, we want to prove the axiom of pairing, that is to say that our construction {a, b} indeed satisfies ∀x (x ∈ {a, b} ⇔ (x ≃ a ∨ x ≃ b)) (lemma 43). In particular, this requires ({a, b}/i(root(a))) ≃ a and ({a, b}/ j(root(b))) ≃ b (lemmas 36 and 37).

We formally prove in our implementation the 53 lemmas. Some of our proofs just follow the informal ones, while some others rely on the type of formulas and its embedding into propositions. The complete proofs can be found at https://github.com/Deducteam/dedukti_set_theory/.

The first two lemmas are immediate consequences of the definition of Leibniz equality in the standard library of DEDUKTI: constant symbol = [ s ] : El s → El s → El prop ; notation = infix 4; rule π ( @ = $s $x $y ) → Π ( P : El $s → El prop ) , π ( P $x ) → π ( P $y );

All the other lemmas of IZmod not involving comp (and therefore restricted formulas) are proved using the blueprint [6, see Proposition 1].

An Example of Proof

To show the way lemmas are proved in DEDUKTI we will take the example of lemma 30 and comment its proof. This lemma states that

a ∈ b ∧ a ≃ c ⇒ c

Lemmas involving Formulas

Lemma 32 is as follows:

(P(z ← a) ∧ a ≃ b) ⇒ P(z ← b)
where P is a restricted formula. Therefore, we formulate lemma 32 using the type of formulas, as well as its interpretation in prop. The valuation update σ z a represents the assignment of variable z ← a. The proof proceeds by induction over formulas: each case is proved easily, using the lemmas that have already been checked by DEDUKTI.

Strong extensionality (lemma 41) is implemented similarly:

Weak Extensionality

The proofs of lemmas 44, 47 and 48 require the use of the axiom of (weak) extensionality [6]. Therefore, we prove in DEDUKTI that extensionality is a consequence of strong extensionality (lemma 41):

Weak extensionality. ∀c∀d (∀z (z ∈ c ⇔ z ∈ d) ⇒ c ≃ d)
Proof. We assume that ∀z (z ∈ c ⇔ z ∈ d). In order to prove c ≃ d, we follow the blueprint given by Dowek and Miquel [5, see Proposition 1] by using an instance of the strong extensionality axiom where R(x, y) is (x ≃ c ∧ y ≃ d) ∨ x ≃ y. We are left to prove the three premises of this instance of strong extensionality:

• (c ≃ c ∧ d ≃ d) ∨ c ≃ d
derives from the reflexivity of bisimilarity (lemma 3).

• We need to prove:

∀a∀a ′ ∀b (a ′ ∈ a ∧ ((a ≃ c ∧ b ≃ d) ∨ a ≃ b) ⇒ (∃b ′ (b ′ ∈ b ∧ ((a ′ ≃ c ∧ b ′ ≃ d) ∨ a ′ ≃ b ′ )) We assume (a ′ ∈ a ∧ ((a ≃ c ∧ b ≃ d) ∨ a ≃ b) and we choose b ′ equal to a ′ , so that (a ′ ≃ c ∧ a ′ ≃ d) ∨ a ′ ≃ a ′

The Axioms of IZst Theory

We have now encoded in DEDUKTI all the axioms of IZst set theory: the strong extensionality axiom corresponds to lemma 41, the axiom of the union is implemented by lemma 42, the pairing axiom corresponds to lemma 43, the axiom of the power set is encoded by lemma 44, the comprehension scheme is implemented by lemma 45, the axiom of infinity corresponds to lemma 51 and the transitive closure axiom is encoded by lemmas 52 and 53.

Conclusion

We have implemented in DEDUKTI a version of set theory -IZst -that corresponds to Zermelo set theory, with the Strong Extensionality axiom and the Transitive Closure axiom. To do so, we have adapted the work by G. Dowek and A. Miquel [START_REF] Dowek | Cut elimination for Zermelo set theory[END_REF] from Deduction modulo theory to λ Π-calculus modulo theory and have encoded sets with a structure of pointed graphs of the IZmod theory. We have formally written all the proofs of the lemmas allowing us to implement set theory in DEDUKTI.

To define and prove the lemmas corresponding to the Comprehension axiom, we have developed a language of formulas, along with operators interpretation and update. In particular, the language of formulas allows us to prove that the Extensionality axiom derives from the Strong Extensionality axiom.

Historically, the encoding of sets by pointed graphs had been designed to enjoy the normalization property. IZmod expressed in Deduction modulo theory does so, but the case of our implementation in λ Π-calculus modulo theory remains to be investigated.

The implementation of IZmod theory represents the first significant corpus of formal proofs in LAMBDAPI. This implementation paves the way for interoperability between DEDUKTI and theorem provers based on set theory, such as ISABELLE/ZF.

  Equality and Membershiprule $a simeq $b → '∃ r : El ( node arrow ( node arrow prop )) , r ( root $a ) ( root $b ) ∧ ( '∀ x , '∀ x ' , '∀ y , eta $a x ' x ∧ r x y ⇒ '∃ y ' , eta $b y ' y ∧ r x ' y ') ∧ ( '∀ y , '∀ y ' , '∀ x , eta $b y ' y ∧ r x y⇒ '∃ x ' , eta $a x ' x ∧ r x ' y '); rule $a ∈ $b → '∃ x , (( eta $b x ( root $b )) ∧ ( $a simeq cr $b x )); Constructions rule eta ( join $a ) $x $x ' → ( '∃ y , '∃ y ' , ( $x = i y ) ∧ ( $x ' = i y ') ∧ eta $a y y ') ∨ ( '∃ y , '∃ z , ( $x = i y ) ∧ ( $x ' = o ) ∧ eta $a y z ∧ eta $a z ( root $a )); rule eta ( pair $a $b ) $x $x ' → ( '∃ y , '∃ y ' , (( $x = i y ) ∧ ( $x ' = i y ') ∧ eta $a y y ')) ∨ ( '∃ y , '∃ y ' , ( $x = j y ) ∧ ( $x ' = j y ') ∧ eta $b y y ') ∨ (( $x = i ( root $a )) ∧ ( $x ' = o )) ∨ (( $x = j ( root $b )) ∧ ( $x ' = o )); rule eta ( powerset $a ) $x $x ' → ( '∃ y , '∃ y ' , ( $x = i y ) ∧ ( $x ' = i y ') ∧ eta $a y y ') ∨ ( '∃ y , '∃ c , ( $x = i y ) ∧ ( $x ' = j (ρ c )) ∧ ( eta $a y ( root $a )) ∧ (( cr $a y ) ∈ c )) ∨ ( '∃ c , ( $x = j (ρ c )) ∧ ( $x ' = o ));rule eta omega $x $x ' → ( '∃ y , '∃ y ' , ( $x = i y ) ∧ ( $x ' = i y ') ∧ ( y < y ')) ∨ ( '∃ y , ( $x = i y ) ∧ ( $x ' = o ) ∧ Nat y ); rule eta ( closure $a ) $x $x ' → ( '∃ y , '∃ y ' , (( $x = i y ) ∧ ( $x ' = i y ') ∧ eta $a y y ')) ∨ ( '∃ y , ( $x = i y ) ∧ ( $x ' = o ) ∧ ( '∀ c : El ( node arrow prop ) , (( '∀ z , eta $a z ( root $a ) ⇒ c z ) ∧ ( '∀ z , '∀ z ' , ( eta $a z z ') ∧ ( c z ') ⇒ ( c z ))) ⇒ c y )); rule root ( join $a ) → o ; rule root ( pair $a $b ) → o ; rule root ( powerset $a ) → o ; rule root omega → o ; rule root ( closure $a ) → o ;

  constant symbol formula : Set ; constant symbol eqF : El nat → El nat → El formula ; constant symbol inF : El nat → El nat → El formula ; constant symbol andF : El formula → El formula → El formula ; constant symbol orF : El formula → El formula → El formula ; constant symbol allF : El nat → El formula → El formula ; constant symbol exF : El nat → El formula → El formula ; constant symbol impF : El formula → El formula → El formula ; constant symbol fF : El formula ; constant symbol tF : El formula ;

  ∈ b Proof. We first assume graphs a, b and c and H the proof of a ∈ b∧a ≃ c. a ∈ b rewrites to ∃x (x η b root b∧ a ≃ (b/x)). We decompose the left part of H as x and Hx which is a proof of x η b root b ∧ a ≃ (b/x). As the goal is to prove c ∈ b, that is, ∃y (y η b root b ∧ c ≃ (b/y)), we need to find a suitable y. We take x and now have two goals: x η b root b and c ≃ (b/x). The first one is proved by applying the left part of Hx. The second one is obtained by applying lemma 5 (transitivity of ≃) to c, a and b/x. To apply lemma 5, we need to prove c ≃ a ∧ a ≃ b/x. c ≃ a is proved by applying lemma 4 (symmetry of ≃) to a, c and the right part of H that is a proof of a ≃ c. a ≃ b/x derives from the right part of Hx. □ This proof is written in DEDUKTI with the following code: opaque symbol lemma30 : π ( '∀ a , '∀ b , '∀ c , (( a ∈ b ) ∧ ( a simeq c )) ⇒ ( c ∈ b )) := begin assume a b c H ; refine ex_e node _ ( and_el _ _ H ) _ _ ; assume x Hx ; refine ex_i node x _ _ ; refine and_i _ _ _ _ { refine and_el _ _ Hx } { refine lemma5 c a ( cr b x ) ( and_i _ _ ( lemma4 a c ( and_er _ _ H )) ( and_er _ _ Hx ))} end ;

  opaque symbol lemma32 : Π ( z : El nat ) , Π ( f : El formula ) , π ( '∀ a , '∀ b , ( '∀ σ : ( El nat → El graph ) , (( interpretation ( update σ z a ) f ) ∧ ( a simeq b )) ⇒ ( interpretation ( update σ z b ) f )))

  opaque symbol lemma41 : Π ( x y : El nat ) , Π ( f : El formula ) , Π ( c d : El graph ) , π ( '∀ σ : ( El nat → El graph ) , (( interpretation ( update ( update σ x c ) y d ) f ) ∧ ( '∀ a , '∀ a ' , '∀ b , (( a ' ∈ a ) ∧ ( interpretation ( update ( update σ x a ) y b ) f )) ⇒ ( '∃ b ' , (( b ' ∈ b ) ∧ ( interpretation ( update ( update σ x a ') y b ') f )))) ∧ ( '∀ b , '∀ b ' , '∀ a , (( b ' ∈ b ) ∧ ( interpretation ( update ( update σ x a ) y b ) f )) ⇒ ( '∃ a ' , (( a ' ∈ a ) ∧ ( interpretation ( update ( update σ x a ') y b ') f ))))) ⇒ ( c simeq d ))

  is a tautology and we are left to prove a ′ ∈ b. If a ≃ b, then by lemma 31 this is a consequence of a ′ ∈ a. Otherwise we have a ≃ c ∧ b ≃ d so with lemma 31 applied twice we are left to prove a′ ∈ c ⇒ a ′ ∈ d, which is a consequence of ∀z (z ∈ c ⇔ z ∈ d).• We proceed similarly for the third premise. □ In order to implement this result in DEDUKTI we first prove the following intermediary lemma:opaque symbol lemmaHypExt : Π ( c d : Graph ) , π (( '∀ z , ( z ∈ c ) ⇔ ( z ∈ d )) ⇒ (((( c simeq c ) ∧ ( d simeq d )) ∨ ( c simeq d )) ∧ ( '∀ a , '∀ a ' , '∀ b , (( a ' ∈ a ) ∧ ((( a simeq c ) ∧ ( b simeq d )) ∨ ( a simeq b ))) ⇒ ( '∃ b ' , (( b ' ∈ b ) ∧ ((( a ' simeq c ) ∧ (b ' simeq d )) ∨ (a ' simeq b '))))) ∧ ( '∀ b , '∀ b ' , '∀ a , (( b ' ∈ b ) ∧ ((( a simeq c ) ∧ ( b simeq d )) ∨ ( a simeq b ))) ⇒ ( '∃ a ' , (( a ' ∈ a ) ∧ ((( a ' simeq c ) ∧ (b ' simeq d )) ∨ (a ' simeq b ')))))))Then in order to prove weak extensionality, we assume graphs c and d, and H the hypothesis ∀z (z ∈ c) ⇔ (z ∈ d). Then we apply lemma 41 to:• natural numbers zero and one• the formula (orF (andF (eqF zero two) (eqF one three)) (eqF zero one))• graphs c and d• the valuation (update (update (λ _, empty_set) two c) three d)• the proof lemmaHypExt c d H.Indeed, in the formula (orF (andF (eqF zero two) (eqF one three)) (eqF zero one)), two will be interpreted by c and three by d, thanks the valuation (update (update (λ _, empty_set) two c) three d). Finally we obtain extensionality: opaque symbol lemmaExt : Π ( c d : El graph ) , π (( '∀ x , ( x ∈ c ) ⇔ ( x ∈ d )) ⇒ ( c simeq d )) := begin assume c d H ; refine lemma41 zero one ( orF ( andF ( eqF zero two ) ( eqF one three )) ( eqF zero one )) c d ( update ( update (λ _ , empty_set ) two c ) three d ) ( lemmaHypExt c d H ) end ;

  El graph ; symbol pair : El graph → El graph → El graph ; symbol powerset : El graph → El graph ; symbol omega : El graph ; // set of natural numbers symbol closure : El graph → El graph ;

symbol eta : El graph → El node → El node → El prop ; symbol root : El graph → El node ; symbol cr : El graph → El node → El graph ; // change of root constant symbol o : El node ; symbol i : El node → El node ; symbol i ' : El node → El node ; symbol j : El node → El node ; symbol j ' : El node → El node ; symbol I : El node → El prop ; symbol J : El node → El prop ; // injection from graphs to nodes ( for powerset ) and its inverse symbol ρ : El graph → El node ; symbol ρ ' : El node → El graph ; // natural numbers constant symbol O : El node ; symbol Succ : El node → El node ; symbol Pred : El node → El node ; symbol Null : El node → El prop ; symbol Nat : El node → El prop ; symbol < : El node → El node → El prop ; symbol simeq : El graph → El graph → El prop ; // bisimilarity symbol ∈ : El graph → El graph → El prop ; symbol join : El graph →

Acknowledgments

The authors would like to thank Frédéric Blanqui for his precious insight on the internal functioning of LAMBDAPI.

Interpretation

In this section, we define the interpretation of the datatype of formulas into the type of propositions. Since formulas can contain variables (represented as natural numbers), the interpretation depends on a valuation from natural numbers to graphs. Hence we have:

The interpretation of quantifiers requires an operator for updating a valuation, taking as arguments the initial valuation, the variable that needs to be updated, and the new value of this variable:

The computation of update σ x a y relies on a decision procedure for the equality on natural numbers: if x = y then update σ x a y is a, other wise it is σ y. Rather than implementing such a decision procedure separately, we embed it in an auxiliary function:

update σ x a y reduces to update1 σ x a y y, and update1 σ x a y z reduces to a if x = y, and σ z otherwise, by decreasing x and y until one (or both) becomes zero. The last argument to update1 is thus used for storing the variable at which update σ x a was initially called. The rewrite rules are as follows:

We can prove that update and update1 satisfy the intended properties:

Now we have all the tools to define the rewriting rules of the interpretation of formulas:

→ '∃ a , interpretation ( update $σ $x a ) $f ; rule interpretation $σ fF → ⊥; rule interpretation $σ tF → ⊤;

Appendix

The complete proofs can be found at https://github.com/Deducteam/dedukti_set_theory/.

Lemma Number of lines in the proof

Lemma Number of lines in the proof


Guy Archambault

Sylvie Gentier

Joëlle Riss

Rock Flamand

Colette Sirieix

A reevaluation of irregular joint shear behavior on the basis of 3D modelling of their morphology -Part Il: Joint Shear Behavior Mechanical Modelling

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Lorick Huang
email: lorick.huang@insa-toulouse.fr.

Mahmoud Khabou
email: mahmoud.khabou@insa-toulouse.fr

The non-linear discrete-time Hawkes process

The non-linear Hawkes process is a point process for which the occurrence of future events depends on the history, either by excitation or inhibition. This property made it popular in many fields, such as neuro-sciences and social-dynamics. In this paper we propose a tractable hidden Markov chain time series as a discrete-time Hawkes process. Our model allows for cross-excitation and inhibition between components, as well as for exogenous random noise on the intensity. We then prove a convergence theorem when the time scale goes to zero. Finally, we suggest a parametric calibration method for the continuous-time Hawkes process based on the regression on the discrete-time approximation.

Introduction

In this paper we introduce a multi-variate time series (P n ) n∈N in N d that represents the bin-count sequence of some events on regular time intervals of size h > 0. Examples of integer valued time series have been studied in the past, such as the g-functions for {0, 1} valued sequences by Berbee [START_REF] Berbee | Chains with infinite connections: uniqueness and Markov representation[END_REF] and the integer valued auto-regressive of order p (INAR(p)) process by Alzaid and Al-Osh [START_REF] Alzaid | An integer-valued pth-order autoregressive structure (INAR(p)) process[END_REF].

The specificity of our time series resides in its hidden Markov chain structure, in which the intensity depends on the history of the process. In fact, for any n ∈ N and any i = 1, • • • , d, the bin-count variable P i n is assumed to be the realisation of a Poisson distribution of parameter hψ i (U n ), where (U n ) n∈N is a predictable tensor Markov chain, i.e. U n ∈ σ(P 0 , • • • , P n-1 , W 0 , • • • , W n-1 ), where (W n ) n∈N are some exogenous random perturbations and (ψ i ) i=1,••• ,d are non-negative functions. The choice of the Poisson distribution for the integer valued random variable is natural: we assume that the time bins are reasonably small, so that one would observe at most one event most of the time, while still allowing for the rare occurrence of two or more events. This, for instance, is not possible with the simpler choice of a Bernoulli distribution (as in [START_REF] Seol | Limit theorems for discrete Hawkes processes[END_REF]), even though that for very small bin size, the two choices are virtually the same. The aforementioned model, called the non-linear discrete-time Hawkes process (by an analogy with the continuous-time Hawkes process made in section 4) allows for self and cross excitation as well as inhibition, which makes it a good candidate for modelling systems in which the occurrence of an event triggers (excitation) or discourages (inhibition) other events.

The main result of this paper is proving that the continuous-time process U ⌊t/h⌋ , ⌊t/h⌋ n=0 P n t∈R+ converges in distribution to a standard continuous-time non-linear Hawkes process, as the time step h goes to zero. The linear uni-variate Hawkes process was introduced in [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] as an example of a self-exciting point process. It was then extended to the more general non-linear multivariate process in [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF]. Initially used in the context of seismic activity in [START_REF] Ogata | Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes[END_REF], Hawkes process found applications in many fields such as genome analysis [START_REF] Reynaud | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF], portfolio credit risk [START_REF] Errais | Affine point processes and portfolio credit risk[END_REF], micro-structure noise [START_REF] Bacry | Modelling microstructure noise with mutually exciting point processes[END_REF] and social networks [START_REF] Bacry | Sparse and low-rank multivariate Hawkes processes[END_REF]. Due to its clear branching structure [START_REF] Hawkes | A cluster process representation of a self-exciting process[END_REF] and the possibility of obtaining closed formulae for its asymptotic moments, the linear Hawkes process has been studied extensively in the literature. Linear Hawkes processes with exponentially decaying intensities have been particularly dealt with in many articles. In fact, the memorylessness of the exponential kernel ensures that the intensity is a Markov process that solves a stochastic differential equation (SDE ). For instance, Giesecke and Kim [START_REF] Giesecke | Estimating tranche spreads by loss process simulation[END_REF] devised a self-exiting point process with CIR intensity (N t ) t∈R+ whose intensity (λ t ) t∈R+ solves the SDE dλ t = β(µ -λ t )dt + σ λ t dW t + αdN t , where (W t ) t∈R+ is a Brownian motion that can be seen as a persistent market noise. Dassios and Zhao [START_REF] Dassios | Efficient simulation of clustering jumps with CIR intensity[END_REF] suggested an efficient simulation procedure for this process, based on the fact that its evolution follows the well-studied CIR dynamics between jumps. Despite the explicit formulae obtained thanks to its Markov properties, the simple exponential kernel is a limiting choice. In fact, it does not account for longer memories nor for delays in the excitation. Similarly, the linearity's tractability comes at an expense: unlike the non-linear Hawkes process, the linear process does not allow for inhibition, nor for bounded excitation. Moreover, the linear process cannot be perturbed by a simple Gaussian diffusion as it can drive the intensity below zero.

The link between time series and linear Hawkes processes have been studied in the seminal works of Kirchner [START_REF] Kirchner | Hawkes and INAR(∞) processes[END_REF], where the linear Hawkes process was proven to be a continuous version of the INAR(∞) time series and vice-versa. The proof of convergence relies on the tightness of the sequence in the space of point processes then showing that all the sub-sequential weak limits follow the same Hawkes dynamics. To the best of our knowledge, the non-linear Hawkes process has no intuitive discrete-time equivalent in the literature. The time series introduced in our paper is a recursively constructed approximation for the Hawkes processes with Erlang kernels (i.e. the product of a polynomial and an exponential) that have a Markov structure. We also allow for the presence of an exogenous Gaussian perturbation that happens continuously in time. We then use specific Markov process techniques to show the weak convergence in the space of càdlàg functions.

The goal of this discrete-time approximation is twofold; on one hand, most data are recorded on regular time intervals (e.g. high-frequency financial data) which makes a discrete-time process more realistic, on the other hand its simulation is more intuitive and tractable than the classic Ogata's thinning algorithm ( [START_REF] Ogata | On Lewis' simulation method for point processes[END_REF]) albeit at the expense of introducing discretisation bias. Moreover, from a practical point of view, the convergence theorem suggests that the parametric calibration of the time-series yields a good estimation of the Hawkes process' parameters, which is verified numerically.

The paper is structured as follows: Section 2 gives the recursive construction of the discrete-time model and explores its Markov structure. In section 3 we give the definition of the non-linear Hawkes process, to which we add a Gaussian perturbation and then we show that it is a Markov process if the kernels are Erlang functions. Section 4 establishes the link between the discrete-time model and the continuous-time process via the convergence of the infinitesimal generators. Finally, in section 5 we run a parametric calibration of the Hawkes process based on a regression on its discrete-time approximation. The experimental results are given with no theoretical guarantees.

2 The discrete-time model

Preliminaries

Throughout this paper, d ∈ N denotes the number of components and h > 0 denotes the width of a time bin. Let ϕ 1 , • • • , ϕ d be a family of non-negative functions, called the jump-rate functions. We set the following assumption

Assumption 2.1. The functions ϕ 1 , • • • , ϕ d are bounded, i.e. ∥ϕ∥ ∞ = max i=1,••• ,d sup s∈R ϕ i (s) < ∞.
For any vector u = (u 1 , • • • , u d ) of non-negative real numbers, we say that the vector ξ follows the multivariate Poisson distribution of parameter u if the variables ξ 1 , • • • , ξ d are independent and

ξ i ∼ Pois(u i ), ∀i = 1, • • • , d.
for a more concise notation, we simply write ξ ∼ Pois(u).

We denote by N q an independent realization of the multivariate standard normal distribution of dimension q.

Construction

For the sake of better understanding, we start this section by building our discrete-time model in the uni-variate setting, with a memory parameter p equal to zero. Knowing the history of the process until n -1, the bin-count sequence P n is supposed to be drawn from a Poisson distribution of parameter hϕ(a n ), where (a n ) n∈N is a the predictable Markov chain that is defined recursively:

a n = µ(1 -e -βh ) + e -βh a n-1 + αP n-1 + σ √ hN (0, 1),
where α ∈ R, σ ≥ 0, β > 0 and ϕ is a positive function, for some fixed intial state. The hidden Markov chain a can be seen as a pre-intensity that takes into account the recent past of the process as well as the exogenous noise (the Gaussian variable), which is turned into the intensity by the action of ϕ.

We simulate the uni-variate bin count sequence P , its cumulative sum H as well as the intensity ϕ((a n )) n∈N (scaled by h) and the pre-intensity (a n ) n∈N for the jump-rate ϕ(x) = ln(1 + e x ) ∧ 40, a decay parameter β = 5, a diffusion parameter σ = 2 and a baseline pre-intensity µ = 2. Figure 1 gives an illustration in the self-exciting case of α = 2 whereas the self-inhibiting case α = -2 is illustrated on figure 2. We now generalize this process to the multivariate setting for a general memory parameter p. This construction can be seen as an Euler scheme of the continuous-time Hawkes process' intensity, which will be defined in Section 3.

Definition 2.2. For an integer p ∈ N, let α = (α ij q ) i,j=1,••• ,d q=0,••• ,p ∈ R (p+1)×d×d . Let β = (β ij ) i,j=1,••• ,d ∈ M d (R * + ), ω ∈ (R * + ) d , σ ∈ R d + and µ ∈ R d . The discrete time Hawkes process H = (H 1 n , • • • , H d n ) n∈N is defined by H n = n k=1 P k ,
where P k ∼ Pois(hϕ(a k )), ∀k ≥ 1 and the "pre-intensities" (a k ) k∈N are defined as

a i k+1 = µ i + d j=1 p q=0 α ij q u ij,q k+1 + σ i y i k+1 ,
where

(u ij,q ) i,j=1,••• ,d q=0,••• ,p
and (y i ) i=1,••• ,d are respectively the discrete cascade of memory terms and the discrete cascade of perturbation terms, defined recursively by:

                       y i k+1 = e -ωih y i k+1 + √ hN (0, 1) u ij,p k+1 = e -βij h (u ij,p k + hu ij,p-1 k ), . . . u ij,q k+1 = e -βij h (u ij,q k + hu ij,q-1 k ), ∀i, j = 1, • • • , d, . . . u ij,0 k+1 = e -βij h u ij,0 k + P j k . (1)
With the given initial states u ij,q 0 = y i 0 = 0, for every i, j = 1, • • • , d and q = 0, • • • , p -1 and P 0 = 0 d . Remark 2.3. For each i = 1, • • • , d the sequence (P n ) can be seen as the number of events of a continuous-time counting process (N t ) t≥0 observed in a time bin of size h, i.e.

P n = N nh -N (n-1)h .
That is why a bin-size parameter h appears in the definition.

The system 1 in Definition 2.2 can be put under the matrix form. Let

U ij k =    u ij,p k . . . u ij,0 k    ∈ R p+1 , Y k =    y 1 k . . . y d k    ∈ R d . (2)
The memory and perturbation cascades are perfectly described by the family

(U ij ) i,j=1,••• ,d , Y . Let α ij = (α ij p , • • • , α ij 0 ) and K =        0 1 • • • 0 0 0 0 • • • 0 0 . . . . . . . . . . . . . . . 0 0 . . . 0 1 0 0 . . . 0 0        ∈ M p+1 (R)
, with 1 on the super-diagonal and 0 elsewhere.

From the definition of K we set M h ij = e -βij h (I + hK), where I is the identity matrix. e p+1 designates the vector of R p+1 with 0 everywhere and 1 in the last component and

Ω h = diag(e -ω1h , • • • , e -ω d h ). Proposition 2.4. The cascade (U n , Y n ) n∈N where U n = (U ij n ) i,j=1,••• ,d is a Markov chain that follows the matrix dynamics U ij n+1 = M h ij U ij n + P j n e p+1 Y n+1 = Ω h Y n + √ hN d Its one step generator T h evaluated at (u, y) ∈ (R p+1 ) d×d × R d for a function f : (R p+1 ) d×d × R d → R has the expression T h f (u, y) =E [f (U n+1 , Y n+1 )| (U n , Y n ) = (u, y)] , =e -h d i=1 ϕi(ã i (u,y)) E N f M h ij u ij i,j=1,••• ,d , Ω h y + √ hN d + e -h d i=1 ϕi(ã i (u,y)) d m=1 E N f M h ij u ij + 1 m=j e p+1 i,j=1,••• ,d , Ω h y + √ hN d hϕ m (ã m (u, y)) + R 2 (u, y).
Where the remainder R 2 is

R 2 (u, y) = E N   i1+•••+i d ≥2 f M h ij u ij + i m e p+1 i,j=1,••• ,d , Ω h y + √ hN d d m=1 (hϕ m (ã m (u, y))) im i m ! e -hϕm(ã m (u,y))   and for every i = 1, • • • , d ãi (u, y) = µ i + σ i g i + d j=1 ⟨α ij , u ij ⟩.
E N stands for the integration with respect to the standard normal distribution in d-dimensions.

Proof. The matrix equation is simply a compact formulation of the recursive definition of (g i , b ij,q ) i,j=1

P i n ∼ P   hϕ i   µ i + σ i g i n + d j=1 p q=0 α ij q b ij,q n     , ∼ P   hϕ i   µ i + σ i g i + d j=1 ⟨α ij , u ij ⟩     , ∼ P hϕ i ãi (u, y) .
The expectation operator E stands for the integration with respect to both the Poisson variable and the Gaussian variable. We denote by E P the sum against the Poisson distributions and E N the expectation with respect to the multivariate standard normal distribution. Since the Poisson variables and the Gaussian variables are built in two separate spaces we have that:

T h f (u, y) =E [f (U n+1 , Y n+1 )| (U n , Y n ) = (u, y)] , =E N E P f (M h ij U ij n + P j n e p+1 ) i,j=1,••• ,d , Ω h Y n + √ hN d |(U n , Y n ) = (u, y) , =E N   +∞ i1,••• ,i d =0 f (M h ij u ij + i j e p+1 ) i,j=1,••• ,d , Ω h y + √ hN d P P 1 n = i 1 , • • • , P d n = i d |(U n , Y n ) = (u, y)   .
Keeping in mind that P i n ∼ P hϕ i ãi (u, y) are independent we have that

T h f (u, y) = +∞ i1,••• ,i d =0 E N   f (M h ij u ij + i j e p+1 ) i,j=1,••• ,d , Ω h y + √ hN d d j=1 hϕ j ãj (u, y) ij i j ! e -hϕj (ã j (u,y))   , =e -h d i=1 ϕi(ã i (u,y)) E N f M h ij u ij i,j=1,••• ,d , Ω h y + √ hN d + e -h d i=1 ϕi(ã i (u,y)) d m=1 E N f M h ij u ij + 1 m=j e p+1 i,j=1,••• ,d , Ω h y + √ hN d hϕ m (ã m (u, y)) + R 2 (u, y).
the last equality is obtained by separating according to the three scenarii

• Exactly zero jumps.

• One jump at exactly one component j ∈ 1, d .

• Two or more jumps.

In order to illustrate the discrete time Hawkes process we simulate it for h = 0.2, d = 2 and p = 1. The jump-rate functions are chosen to be ϕ 1 (x) = max(x, 0) ∧ 40 and ϕ 2 (x) = log(1 + e x ) ∧ 40. The decay rates are chosen to be β 11 = β 12 = 5, β 21 = β 22 = 6 and the baseline parameters are (µ 1 , µ 2 ) = (2, 3). The results are shown on figures 3 and 4 in the appendix 7.

Notice how H 1 (or H 2 ) has jumps whose sizes are larger than one. This is possible because the jumps are chosen to be Poisson variables (along the lines of [START_REF] Kirchner | An estimation procedure for the Hawkes process[END_REF]), unlike the alternative choice of a Bernoulli variable (cf. [START_REF] Seol | Limit theorems for discrete Hawkes processes[END_REF]). Even though the two choices are equivalent when h → 0, the Poisson approach allows multiple jumps within a single time bin and allows the possibility hϕ i (a i k ) > 1.

Remark 2.5. Even though in Proposition 2.4 the cascade (U n , Y n ) n∈N were proven to be a Markov chain, it is straightforward to show that (U n , Y n , H n ) n∈N is also a Markov chain. From now on, we mainly focus on the pre-intensity, because it already contains the jumps of H.

Stability

We show that if the jump rates are bounded, the pre-intensities

U ij k i,j=1,••• ,d k∈N
have bounded averages and are thus almost surely finite.

Proposition 2.6. Let (U n , Y n ) n∈N be the Markov chain defined in 2 following the recursive formula of Proposition 2.4 and assume assumption 2.1 is in force. Then there is a positive constant C independent of n such that

E   d i,j=1 p q=0 u ij,q n + |y i n |   ≤ C
for any n ∈ N.

Proof. We start with Y n , using the triangular inequality and the mean of the folded normal distribution

E[|y i k+1 |] ≤ e -ωih E[|y i k |] + 2h π ,
hence, using the discrete version of Gronwall's lemma:

E[|y i n+1 |] ≤ 2h π 1 1 -e -ωih .
For U n , we set the decreasing positive sequence (π ij 0 , • • • , π ij q ) and we define

D ij k+1 = p q=0 E[u ij,q k+1 ]π ij q β q+1 ij . The dynamics of U ij k+1 yield D ij k+1 = e -βij h p q=0 E u ij,q k π ij q β q+1 ij + h p q=1 E u ij,q-1 k π ij q β q+1 ij + E P j k π ij 0 β ij , = e -βij h E u ij,p k π ij p β p+1 ij + p-1 q=0 E u ij,q k π ij q + hβ ij π ij q+1 β q+1 ij + E P j k π ij 0 β ij , ≤ e -βij h E u ij,p k π ij p β p+1 ij + p-1 q=0 E u ij,q k π ij q + hβ ij π ij q+1 β q+1 ij + h∥ϕ∥ ∞ π ij 0 β ij .
Using the fact that 1 + x ≤ e x we get

e -βij h π ij q + hβ ij π ij q+1 = π ij q e -βij h 1 + hβ ij π ij q+1 π ij q , ≤ π ij q exp hβ ij π ij q+1 π ij q -1 ,
where

π ij q+1 π ij q -1 < 0.
Finally, we obtain that

D ij k+1 ≤ κ ij D ij k + h∥ϕ∥ ∞ π ij 0 β ij with κ ij = max q=0,••• ,p-1 exp hβ ij π ij q+1 π ij q -1 , e -βij h < 1,
hence the result.

The result above does not prove stability in variation (cf. [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF]), but it ensures that the self-excitation (in case α is positive) is not explosive. In fact, the result of Proposition 2.6 entails that (using Markov's inequality)

P u ij,q n ≥ M ≤ C M , ∀n ∈ N.
The result also ensures that the sequence (U n , Y n ) n∈N is tight, thus having weakly convergent sub-sequence, which in turn can be used to prove the convergence towards a stationary distribution. We chose not to prove it in this paper, as the long-time behaviour has already been studied in [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF] for the continuous-time Hawkes process, using classical Markov process techniques.

3 The continuous-time Hawkes process

Preliminaries

Let Ω d be the space of configurations, where

Ω := ω = n i=1 δ (ti,θi) , 0 = t 0 < t 1 < • • • < t n , θ i ∈ R + , n ∈ N ∪ {+∞} .
Let F be the σ-field associated to the vague topology on Ω d , and P Γ the Poisson measure under which the family

Γ = Γ j j=1,••• ,d where Γ j ([0, t] × [0, b]) (ω j ) := ω j ([0, t] × [0, b]) , (t, b) ∈ R 2 + , j = 1, • • • , d is
a family of independent homogeneous Poisson processes with intensity measures dt⊗dθ. We set

F Γ = F Γ t t≥0
to be the natural filtration of Γ j j=1,••• ,d . The expected value with respect to the Poisson measure is denoted by E Γ . On a separate space Ω ′ , we construct a standard d-variate Brownian motion W = (W 1 , • • • , W d ) and we denote its filtration by F W = F W t t≥0 , and by P W the associated probability. The expected value with respect to the Brownian measure is denoted by E W . We define P = P Γ ⊗ P W and F = F Γ ∨ F W . The expected value with respect to P is denoted by E.

Definition and construction

We start by defining the non-linear multivariate Hawkes process in the general framework.

For i = 1, • • • , d, consider the events times τ i 1 , τ i 2 ,
• • • associated with the i-th component and define the counting process

N i t = i≥1 1 τ i i ≤t . To N = (N i ) i=1,••• ,d we associate a predictable intensity vector λ = (λ i ) i=1,••• ,d such that P N i t+dt -N i t = 1|F t-= λ j t dt
which measures how likely it is for N j to jump between t and t + dt, right before t. Let (δ ij ) i,j=1,••• ,d be a set of real integrable (not necessarily non-negative) functions on R + and ϕ as defined in subsection 2. The process N is called a non-linear Hawkes process of memory kernel δ and jump rate ϕ if its intensity follows the dynamics

λ i t = ϕ i   µ i + [0,t)   d j=1 δ ij (t -s)dN j s   + σ i e -ωi(t-s) dW i s   , = ϕ i   µ i + d j=1 τ j n <t δ ij (t -τ j n ) + [0,t) σ i e -ωi(t-s) dW i s   ,
for some ω i > 0, for every i = 1, • • • , d. [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF] constructed the non-linear Hawkes process as the result of the embedding from a Poisson measure. We give the result in the following theorem.

∈ R d , ω ∈ (R * + ) d and σ ∈ R d + .
The following SDE has a unique solution (N , λ) with N F-measurable and λ F-predictable

         N i t = (0,t]×R+ 1 {θ≤λ i s } Γ i (ds, dθ), t ≥ 0, i = 1, • • • , d, λ i t = ϕ i µ i + [0,t) d j=1 δ ij (t -s)dN j s + σ i e -ωi(t-s) dW i s t ≥ 0 i = 1, • • • , d. (3)
Remark 3.2. In [START_REF] Brémaud | Stability of nonlinear Hawkes processes[END_REF], the assumption that the functions ϕ 1 , • • • , ϕ d are bounded is relaxed and replaced by the condition that ϕ i is l i -Lipschitz and that the matrix

l i +∞ 0 |δ ij (s)|ds i,j=1,••• ,d
has a spectral radius strictly less than one. This guarantees that N exists on R + and has a stationary distribution.

In our definition a Gaussian perturbation was added to the original definition of the intensity. This can be seen as the effect of exogenous factors on each component of the system, that do not occur on the jump times of a point process, but continuously in time. On the other hand, the unperturbed Hawkes intensity follows a piece-wise deterministic dynamics. In the following section, we choose the memory kernels from a specific family of functions and we show how this choice guarantees that the Hawkes process is Markov, without being too restrictive.

Markov properties of the Hawkes process with memory kernels given by an Erlang function

From now on, we restrict the choice for the memory kernels to be in the Erlang functions family

δ ij (s) = p q=0 α ij q s q q! e -βij s ,
for all s ≥ 0 and ∀i, j = 1, • • • , d. Thus, the intensity takes the form

λ i t = ϕ i   µ i + d j=1 [0,t) e -βij (t-s) p q=0 α ij q (t -s) q q! dN j s + σ i [0,t) e -ωi(t-s) dW i s   .
When the memory kernels are Erlang functions, the stability condition in remark 3.

l i p q=0 |α ij 0 | β q+1 ij i,j=1,••• ,d
has a spectral radius strictly less than one.

Remark 3.4. This is a special class of Erlang memory kernels. The more general approach would be to make the (β i ) i=1,••• ,d dependent on q, along the lines of [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF]. This choice would allow to approximate every L 1 (R + , R) kernel, hence approximating any Hawkes process with an integrable kernel. In fact, it has been proven in [START_REF] Mads Bonde Raad | Stability and mean-field limits of age dependent Hawkes processes[END_REF] that for any T there is a constant C T such that for two Hawkes processes N and Ñ of memory kernels δ and δ respectively, we have that

d i=1 E ∥N i -Ñ i ∥ T V ≤ C T d i,j=1 T 0 |δ ij (s) -δij (s)|ds, where ∥ • ∥ T V is the total variation distance on [0, T ].
The results in this paper are entirely adaptable to those kernels, we chose however to write them down solely for β independent from q because the computations are less tedious.

We now show that, the vector formed by the continuous-time memory and perturbation cascades

v ij,q t := [0,t) (t -s) q q! e -βij (t-s) dN j s , (4)
and

z i t := [0,t) e -ωi(t-s) dW i s , (5)
is a (p + 1)d 2 + d-variate Markov process that solves a certain SDE.

Theorem 3.5. Let N be a d-variate Hawkes process with Erlang memory kernels. For every i = 1,

• • • , d, set α ij = (α ij p , • • • , α ij 0 ) ∈ R p+1 . Let V ij t = (v ij,p t , • • • , v ij,0 t ) ∈ R p+1 and Z t = (z 1 t , • • • , z d t ) ∈ R d , such that λ i t = ϕ i   µ i + σ i z i t + d j=1 ⟨α ij , V ij t ⟩   , = ϕ i ãi (V t , Z t ) . Then, (V , Z) = (V ij ) i,j=1,••• ,d , (z i ) i=1,••• ,d is a
Markov process that follows the SDE:

dV ij t = (K -β ij I)V ij t dt + dN j t e p+1 , dZ t = -diag(ω i , • • • , ω d )Z t dt + dW t
where e p+1 = (0,

• • • , 0, 1) and W = (W 1 , • • • , W d ).
Proof. For the auxiliary processes v ij,k for fixed i, j ∈ 1,

• • • , d and q ∈ 1, • • • , p we have v ij,q t e βij t = 1 q! t 0 (t -s) q e βij s dN j s , = 1 q! q l=0
C l q t l t 0 (-s) q-l e βij s dN j s , which yields after differentiation

dv ij,q t + β ij v ij,q t dt e βij t = 1 q! dt q l=0 C l q lt l-1 t 0 (-s) q-l e βij s dN j s + q l=0
C l q t l t(-t) q-l e βij t dN j t , and since

C l q l = q! (q -l)!(l -1)! = q (q -1)! (q -1 -(l -1))!(l -1)! = qC l-1 q-1 ,
we have

dv ij,q t + β ij v ij,q t dt e βij t = 1 q! dt q l=1 qC l-1 q-1 t l-1 t 0 (-s) q-l e βij s dN j s , = 1 (q -1)! dt q-1 l=0 C l q-1 t l t 0 (-s) q-1-l e βij s dN j s , = 1 (q -1)! dt t 0 (t -s) q-1 e βij s dN j s ,
and after re-arranging the terms

dv ij,q t = -β ij v ij,q t + t 0 (t -s) q-1 e -βij (t-s) dN j s dt, = -β ij v ij,q t + v ij,q-1 t dt.
For the last term

v ij,0 t = t 0 e -βij (t-s) dN j t , a straightforward differentiation yields dv ij,0 t = -β ij v ij,0 t dt + dN j t .
The SDE for D i is obtained identically, in fact, it is a multi-variate Ornstein-Uhlenbeck process.

Remark 3.6. Note that the vector process formed by the memory and perturbation cascades as well as the point process (V , Z, N ) is also a Markov process. But since it makes the notations and the computations more tedious, we simply consider (V , Z) that already contains all the jump information. In fact, the jumps of N i are exactly the points of the underlying Poisson measure under the curve ϕ i ãi (V t , Z t ) :

dN i t = R+ 1 {θ≤ϕi(ã i (V t,Zt))} Γ i (dt, dθ).
This remark is the continuous-time analogue of Remark 2.5.

In the next proposition we give the infinitesimal generator for the process (V , Z) defined in the last theorem.

Proposition 3.7. We recall that e p+1 is the vector that has 1 at the p + 1-th entry and 0 elsewhere. Let

f : (R p+1 ) d×d × R d → R (v, z) → f (v, z).
be a C 1 function in the first variable and C 2 in the second variable. Then the infinitesimal generator A of the process (V , Z) at f is

Af (v, z) = d i=1 ∂ 2 f ∂z 2 i (v, z) + d i,j=1 (K -β ij I)v ij , ∇ vij f (v, z) -⟨diag(ω 1 , • • • , ω d )z, ∇ z f (v, z)⟩ + d m=1 ϕ m (ã m (v, z)) (f ((v ij + 1 j=m e p+1 ) i,j=1,••• ,d , z) -f (v, z))
,

where v = (v ij ) i,j=1,••• ,d , z = (z 1 , • • • , z d ) and ãm (v, z) = µ m + σ m z m + d j=1 ⟨α mj , v mj ⟩. Proof.
For the unperturbed uni-variate case, the expression of the generator can be found in [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF]. The second derivatives are a result of the Gaussian perturbation. This is more precisely the generator of a multi-variate Ornstein Uhlenbeck process.

For the linear Hawkes process, (ϕ i (x) = x, ∀i = 1, • • • , d and δ ij ≥ 0, ∀i, j = 1, • • • , d) the Gaussian perturbation with constant volatilities poses a problem, since the intensity is no longer guaranteed to remain non-negative. That is why the volatility must be a variable function in the intensity. For instance, the intensity of an uni-variate linear Hawkes process with an exponential kernel δ(t) = αe -βt is assumed to solve the SDE dλ t = β(µ -λ t )dt + αdN t + σ λ t dW t to ensure that he intensity remains positive (cf. [START_REF] Giesecke | Estimating tranche spreads by loss process simulation[END_REF] and [START_REF] Dassios | Efficient simulation of clustering jumps with CIR intensity[END_REF]), with the upper bound σ 2 ≤ 2µβ. The existence of these "jump-diffusion" processes as well as other results on the affine point processes can be found in [START_REF] Duffie | Transform analysis and asset pricing for affine jumpdiffusions[END_REF].

Convergence of the discrete-time Hawkes process to the continuoustime Hawkes process 4.1 Preliminaries

Let D(R + , E) be the space of right continuous processes with left limits (càdlàg) on a metric space E. We denote by " ⇒ " the weak convergence of processes in the Skorokhod topology (cf. [START_REF] Billingsley | Convergence of probability measures[END_REF] for example). Throughout this paper, Ĉ(E) denotes the set of continuous functions from E to R that vanish at infinity and the functional convergence is in the uniform norm denoted by ∥ • ∥. In the Markov framework, the weak convergence of a family of processes in the Skorokhod topology can be deduced directly from the behavior of its transition operators. This is formulated in the following theorem in [START_REF] Stewart | Markov processes[END_REF]: Theorem 4.1. Let E be locally compact and separable. For h > 0, let ν h (x, •) be a transition function on E × B(E) such that the operator T h defined by

T h f (x) = E f (y)ν h (x, dy)
satisfies T h : Ĉ(E) → Ĉ(E). Suppose that T (t) is a Feller semi-group on Ĉ(E).

Assume that h ↓ 0 and suppose that for every f ∈ Ĉ(E),

lim h↓0 T ⌊t/h⌋ h f = T (t)f, t ≥ 0.
for each h > 0, let Y h k k∈N be a Markov chain on E with a transition function ν h , and suppose that (Y h 0 ) h>0 has a limiting distribution when h ↓ 0. Define

X h t = Y h ⌊t/h⌋ .
Then there is a Markov process X corresponding to T (t) with sample paths in D(R + , E) and

X h ⇒ X.
The the convergence for the semi-groups is usually proved by proving the convergence of the infinitesimal generators on a subset of functions called "the core" of the limiting generator. For a generator A of a contraction semi-group T (t) defined on a Banach functional space L, we suppose that there are two dense subsets D and D 0 such that D 0 ⊂ D ⊂ L. We say that D is a core for A if for any t ≥ 0, T (t) : D 0 → D. Proposition 4.2. Let L be a functional Banach space. For h > 0, let T h be a linear contraction on L and put

A h = T h -I h ,
where I stands for the identity operator. Let T (t) be a strongly continuous contraction semi-group on L with generator A. Let D be a core of A, then the following are equivalent

1. For each f ∈ L, lim h↓0 T ⌊t/h⌋ h f = T (t)f, ∀t ≥ 0. 2. For each f ∈ D, there exists a sequence f h ∈ L such that lim h↓0 f h = f and lim h↓0 A h f h = Af .
Proof. This is a modification of theorem 6.5 in [START_REF] Stewart | Markov processes[END_REF] page 31.

Remark 4.3. Even though the auxiliary processes defined in 4 are predictable (hence left continuous with right limits or càglàd), the theorems used in this section prove the convergence in D(R + , E). That is why we replace our processes with their càdlàg version that has the exact same jumps and dynamics.

Main result

In this section we prove that the discrete time Hawkes process defined in section 2.2 is an approximation of the continuous-time multivariate non-linear Hawkes process. To do so we provide a convergence theorem for the auxiliary processes that contain the jumps and the dynamics. Note that it is possible to prove the results for the augmented process (V , Z, N ), but this is not necessary as (V , Z) is sufficient to reconstruct N (cf. Remark 3.6).

For the initial state, there is nothing to prove since we take it always to be equal to zero. To prove the convergence of the process on D(R + , E) we state the following result Proposition 4.4. Let T h be the discrete cascade's transition operator defined in 2.4 and let A be the generator for the auxiliary processes of a Hawkes process defined in 3.7. Set

A h = T h -I h .
Then, for any f ∈ S(E) (the Schwartz space of functions) one has

lim h↓0 A h f = Af.
Proof. Let f ∈ S(E) and h > 0. Using Lemma 7.2 we have that

T h (u, y) =e -h d i=1 ϕi(ã i (u,y)) E N f M h ij u ij i,j=1,••• ,d , Ω h y + √ hN d + e -h d i=1 ϕi(ã i (u,y)) d m=1 E N f M h ij u ij + 1 m=j e p+1 i,j=1,••• ,d , Ω h y + √ hN d hϕ m (ã m (u, y)) + O(h 2 ).
Note that the subscript i, j = 1, • • • , d has been omitted to make the notations less tedious. For a d-variate Gaussian variable N d , we write a Taylor expansion in the second component

f M h ij u ij , Ω h y + √ hN d =f M h ij u ij , Ω h y + √ h N d , (∇ y f ) M h ij u ij , Ω h y + h 2 N d , (Hess y f ) M h ij u ij , Ω h y N d + h 3/2 |r|=3 R r M h ij u ij , Ω h y (N d ) r ,
where the remainder R r satisfies

|R r M h ij u ij , Ω h y | ≤ max |s|=3 sup z∈R d |D s f (z)|, ≤ C,
since f is in the Schwartz space. Hence by taking the expected value with respect to the centered d-variate Gaussian we have that

E N f M h ij u ij , Ω h y + √ hN d = f M h ij u ij , Ω h y + h 2 d i=1 ∂ 2 f ∂y 2 i M h ij u ij , Ω h y + O(h 3/2 ).
Now, a second order Taylor expansion with a Lagrange remainder in the first component yields

f M h ij u ij , Ω h y =f u ij , Ω h y + d i,j=1 (M h ij -I)u ij , ∇ uij f (u ij , Ω h y) + d i,j=1 |s|=2 1 s! D s f ((1 -c)I + cM h ij )u ij , Ω h y (M h ij -I)u ij s
for some c ∈ (0, 1). Since for any i, j = 1,

• • • , d M h ij -I = h(K -β ij I) + O(h 2 ), f M h ij u ij , Ω h y =f u ij , Ω h y + h d i,j=1 (K -β ij I)u ij , ∇ uij f (u ij , Ω h y) + d i,j=1 O(h 2 )u ij , ∇ uij f (u ij , Ω h y) + h 2 d i,j=1 |s|=2 1 s! D s f ((1 -c)I + cM h ij )u ij , Ω h y ((K -β ij I + O(h))u ij ) s .
Since f is taken to be in the Schwartz space, it follows that its derivatives multiplied by any polynomial are uniformly bounded, hence

f M h ij u ij , Ω h y = f u ij , Ω h y + h d i,j=1 (K -β ij I)u ij , ∇ uij f ((u ij ), Ω h y) + O(h 2 )
A similar argument for the second component of the function and its gradient gives us that

f M h ij u ij , Ω h y =f ((u ij ), y) + h d i,j=1 (K -β ij I)u ij , ∇ uij f ((u ij ), y) -h ⟨diag(ω 1 , • • • , ω d )y, ∇ y f ((u ij ), y)⟩ + O(h 2 ).
Applying the same arguments to the second derivative yields

E N f M h ij u ij , Ω h y + √ hN d =f ((u ij ), y) + h d i,j=1 (K -β ij I)u ij , ∇ uij f ((u ij ), y) -h ⟨diag(ω 1 , • • • , ω d )y, ∇ y f ((u ij ), y)⟩ + d i=1 ∂ 2 f ∂y 2 i u ij , y + O(h 3/2 ).
Finally, using the fact that the jump rates ϕ 1 , • • • , ϕ d are bounded and the Taylor expansion of the exponential

e -h d i=1 ϕi(ã i (u,y)) = 1 -h d i=1 ϕ i ãi (u, y) + O(h 2 ) we conclude that T h (t)f ((u ij ), y) = f ((u ij ), y) + hAf ((u ij ), y) + O(h 3/2 ),
which yields the result.

We now state the main result, establishing the convergence of the discrete-time Hawkes process to the continuous-time Hawkes process. The convergence is illustrated in the transition from figure 5 (h = 0.125) to figure 6 (h = 5.10 -3 ).

Parameters calibration using regression

The estimation of linear Hawkes processes have been thoroughly studied in the linear case, both in a parametric and a non-parametric way.

For the exponential kernels, [START_REF] Ozaki | Maximum likelihood estimation of Hawkes' self-exciting point processes[END_REF] proved that the computation likelihood function is linear in the number of events, numerical results are given in the 1-d case as well. [START_REF] Da | Hawkes Process: Fast Calibration, Application to Trade Clustering, and Diffusive Limit[END_REF] derived an estimation of the parameters using the stationary moments that can be explicitly derived for the linear Hawkes process with exponential kernels. This method is not applicable in the non-linear framework since the expressions of the first moments are not known explicitly. Since the kernel's shape is usually not known in advance and since the exponential family is not rich enough (e.g. power laws, non monotonous functions), methods for non-parametric estimation have been derived for the linear Hawkes process. The first of which is based on the numerical resolution of a system of Wiener-Hopf equations, introduced by [START_REF] Bacry | First-and second-order statistics characterization of Hawkes processes and non-parametric estimation[END_REF]. The second is introduced by [START_REF] Kirchner | An estimation procedure for the Hawkes process[END_REF] and by [START_REF] Eichler | Graphical modeling for multivariate Hawkes processes with nonparametric link functions[END_REF] based on time series. These methods are close to the ones we present in this paper. The reader can refer to [START_REF] Shlomovich | A Parameter Estimation Method for Multivariate Aggregated Hawkes Processes[END_REF] for a comparison between the time series method and an E-M method.

The non-linear framework has also been studied, though to a lesser extent. For the parametric estimation, [START_REF] Lemonnier | Nonparametric Markovian Learning of Triggering Kernels for Mutually Exciting and Mutually Inhibiting Multivariate Hawkes Processes[END_REF] derived an algorithm for the log-likelihood maximisation for kernels composed by an exponential mixture and a ReLU jumps rate. For a more general choice of the jump rate functions as well as kernels, [START_REF] Sulem | Bayesian estimation of nonlinear Hawkes process[END_REF] obtained concentration rates of the posterior distribution on the parameters. When it comes to the non parametric approach, Reynaud-Bouret et al. devised a method for the non-parametric estimation for multivariate Hawkes processes with ReLU jump-rates.

In this section, we assume that both assumptions 2.1 (boundedness of the jump-rate) and 3.3 (jump-rate functions are l i -Lipschitz and the matrix l i p q=0

|α ij 0 | β q+1 ij i,j=1,••• ,d
has a spectral radius less than one) are in force. In fact, if assumption 3.3 is not respected, the pre-intensity ãi can grow to infinity, before being capped by the upper bound ∥ϕ i ∥ ∞ , transforming the Hawkes process into a Poisson process of intensity ∥ϕ i ∥ ∞ from which no other parameters can be recovered. This is the problem of identifiability which is discussed in more details in [START_REF] Sulem | Bayesian estimation of nonlinear Hawkes process[END_REF]. We then implement a conditional least-squares (CLS) regression on the observed discrete time version

(P k ) k=1,••• ,M = (H k -H k-1 ) k=1,••• ,M of an underlying Hawkes process (N t ) t∈[0,T ] , with time bins h = T /M . Let F k = σ (P 1 , • • • , P k )
be the filtration that contains the information up to k. We seek to find min

(α,β,µ)∈(R (p+1)×d×d )×M d (R+)×R d + M k=1 ∥P k -E [P k |F k-1 ]∥ 2 .
For the sake of brevity, we set

Q M = M k=1 ∥P k -E [P k |F k-1 ]∥ 2 , = M k=1 d i=1 P i k -E P i k |F k-1 2 .
The underlying process is supposed to follow multivariate Hawkes dynamics

λ i t = ϕ i   µ i + d j=1 [0,t) δ ij (t -s)dN j s + σ i [0,t) e -ωi(t-s) dW i s   ,
with known jump-rates ϕ, unknown baseline intensities µ and unknown kernels δ ij (s) = e -βij p q=0 α ij q s q q! , where p is unknown. We do not estimate the Gaussian perturbation's parameters. Since we have proved that the discrete time Hawkes process converges to the continuous-time Hawkes process, we set

P i k = N i hk -N i h(k-1)
for a fixed h > 0.

General results about CLS estimators for time-series, including asymptotic bias and normality can be found in the seminal work of [START_REF] Klimko | On conditional least squares estimation for stochastic processes[END_REF].

Using the definition of the discrete time Hawkes process, Q M can be computed linearly in the number of points as

Q M +1 = Q M + d i=1 P i M +1 -E P i M +1 |F M 2 , = Q M + d i=1   P i M +1 -hϕ i   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     2 ,
where the U ij are also computed recursively

U ij M +1 = M ij U ij M + P j M e p+1 .
If the jump-rates are assumed to be differentiable almost everywhere, the same holds for the gradients which are computed recursively

∂ µi Q M +1 = ∂ µi Q M -2h(ϕ i ) ′   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     P i M +1 -hϕ i   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     , ∇ α ij Q M +1 = ∇ α ij Q M -2hU ij M (ϕ i ) ′   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     P i M +1 -hϕ i   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     , ∂ βij Q M +1 = ∂ βij Q M -2h α ij , ∂ βij U ij M +1 (ϕ i ) ′   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     P i M +1 -hϕ i   µ i + d j=1 ⟨α ij , U ij M +1 ⟩     , ∂ βij U ij M +1 = M ij ∂ βij U ij M -hU ij M .
Note that unlike the linear Hawkes process with Erlang kernels, this cannot be simplified to a computation in the number of events. As a first trial we simulate an unperturbed 3-dimensional DTHP with h = 0.1 for M = 5000. The kernels are

δ(s) =  
(1 + 0.2s)e -3s (0.2 + 3s)e -4s 3se -3s 0 (-0.1 -2s)e -2s

(1 -5s)e -1.5s (-0.5 -2s)e -3s 0 (0.7 + 6.2s)e -3.5s The minimization of Q M is done using the SLSQP method of the function scipy.optimize.minimize with assumption 3.3 as a constraint. The initial guesses are αij q ∼ U([0, 1]), βij = 1 and μi = ϕ i -1

M k=1 P i k T .
The minimization results are then averaged (for five calibrations) and are given in figure 7.

Calibration experiment for the perturbed Hawkes process

The underlying dynamics are assumed to follow a continuous-time 3-variate Hawkes process with a Gaussian perturbation of parameters ω = 3 4 5 and σ = 1 2 2 .

We simulate 35 Monte Carlo samples for t ∈ [0, 100], which we regroup in 35 DTHP arrays with a time step h = 0.25. The discretisation introduces bias, by discarding the possibility that a given event in a given bin be the origin of a later event in the same bin. The result of the CLS minimisation is given in figures 8, 9 and 10.

Calibration with a lower degree

We now simulate 15 continuous-time Hawkes processes with kernels δ(s) =   (0.7 -2s + 1.5s 2 )e -1.5s (0.5 + 4s -10s 2 )e -3s 2s 2 e -2s 0 (5s + 6s 2 )e -3s (1 -10s 2 )e -2.5s -5s 2 e -2s 0 (0.8 + 7s + 8s 2 )e -5s The processes are then discretized on a grid with a time-step h = 0.25. We minimise the conditional least squares Q M with assumed kernels of the form

δ ij (s) = α ij 0 e -βij s .
The exponential kernels composed with the average parameters are illustrated in figure 11. When it comes to the baseline intensities refer to figure 12 6 Conclusion

In this paper, we introduced a tractable and intuitive Markov Chain approximation of the multivariate nonlinear Hawkes process with a Gaussian perturbation, and we proved that this approximation converges when the time-bin size goes to zero. This result can be extended to a wide array of similar processes with Markov dynamics, for instance for a perturbation parameter σ that depends on the system's state, a baseline intensity formed by a mixture of Erlang and trigonometric functions to account for seasonality or for jumps with random sizes.

A parametric calibration method has been studied numerically but without theoretical guarantees on the convergence of the estimated parameters to the ground truth parameters. The experiments show that the calibrated parameters are close to the actual kernels despite the Gaussian perturbations and the bias introduced by the discretization. This method can be enhanced for instance by adding a penalisation or by decreasing the variance.

Lemmata

Lemma 7.1. Set E = (R p+1 ) d×d × R d . Then the space Ĉ(E) of continuous functions that vanish at infinity equipped with the uniform norm is a Banach space. Moreover, the space of C ∞ functions of compact support C ∞ c (E) is dense in Ĉ(E). Proof. Let (f n ) n∈N be a Cauchy sequence in Ĉ(E). For every x ∈ E, the sequence (f n (x)) n∈N satisfies

|f n (x) -f m (x)| ≤ ∥f n -f m ∥, ∀n, m ∈ N,
and thus, is a Cauchy sequence in E. Since E is complete, the sequence converges to the point-wise limit f (x). Let p ∈ N and x ∈ E

|f p (x) -f (x)| = |f p (x) -lim n→+∞ f n (x)|, = lim n→+∞ |f p (x) -f n (x)|, ≤ lim n→+∞ ∥f p -f n ∥. Let ϵ > 0 and M such that for m, n ≥ M , ∥f m -f n ∥ ≤ ϵ. If p ≥ M |f p (x) -f (x)| ≤ ϵ.
Thus (f p ) p∈N converges to f uniformly, hence f is continuous. For some p such that ∥f p -f ∥ ≤ ϵ/2, Let K p be the compact such that ∀x / ∈ K p :

|f p (x)| ≤ ϵ/2.
We have for any

x / ∈ K p |f (x)| ≤ |f (x) -f p (x)| + |f p (x)|, ≤ ∥f -f p ∥ + ϵ/2, ≤ ϵ.
Which proves that Ĉ(E) is a Banach space.

For the density of the C ∞ c (E), it is an application of Stone-Weierstrass' theorem, thanks to the existence of bump functions. Lemma 7.2. Assume that assumption 2.1 is in force. Let R 2 be the remainder defined in Proposition 2.4. Then, uniformly in u and y,

R 2 (u, y) = O(h 2 ). Proof. We recall that for a function f ∈ D, R 2 (u, y) = E N   i1+•••+i d ≥2 f M h ij u ij + i m e p+1 i,j=1,••• ,d , Ω h y + √ hN d d m=1 (hϕ m (ã m (u, y))) im i m ! e -hϕm(ã m (u,y))   .
Since f is bounded we have that

|R 2 (u, y)| ≤ E N   i1+•••+i d ≥2 ∥f ∥ ∞ d m=1 (hϕ m (ã m (u, y))) im i m ! e -hϕm(ã m (u,y))   , = ∥f ∥ ∞ i1+•••+i d ≥2 d m=1 (hϕ m (ã m (u, y))) im i m ! e -hϕm(ã m (u,y)) . Since i 1 + • • • + i m ≤ 2 and h < 1, h im ≤ h 2 for any m = 1, • • • , d. Hence |R 2 (u, y)| ≤ ∥f ∥ ∞ h 2 i1+•••+i d ≥2 d m=1 ϕ m (ã m (u, y)) im i m ! e -hϕm(ã m (u,y)) , ≤ ∥f ∥ ∞ h 2 i1+•••+i d ≥2 d m=1 ∥ϕ m ∥ im i m ! , ≤ ∥f ∥ ∞ h 2 e d m=1 ∥ϕm∥ . Lemma 7.3. For i, j = 1, • • • , d, let M h ij = e -βij h (I + hK) and Ω h = diag(e -ω1h , • • • , e -ω d h ) be the matrices defined in section 2.2. Then M h ij = I + h(K -β ij I) + O(h 2 ), and
Ω h = I -diag(ω 1 , • • • , ω d )h + O(h 2 ).
Proof. This is a simple second order Taylor expansion.

Lemma 7.4. Set E = (R p+1 ) d×d × R d .
Let T be the semi-group associated with the Markov process defined by the auxiliary processes 4 and 5, i.e. defined for every t ≥ 0 and f ∈ Ĉ(E) by

T (t)f (v, z) = E [f (V t , Z t )|(V 0 , Z 0 ) = (v, z)] .
Then T is a Feller semi-group.

Proof. It is trivial to see that for any t ≥ 0, T (t) is a contraction (i.e. ∥T (t)f ∥ ≤ ∥f ∥ for every f ∈ Ĉ(E) ) and positive (i.e. for every non-negative function f , T (t)f is a non-negative function). Now we show that it is strongly continuous. To do so, one must show that it has the Feller property

∀f ∈ Ĉ(E), t ≥ 0, T (t)f ∈ Ĉ(E).
First, by solving the SDE in 3.5 with the initial conditions V ij 0 = v ij and Z 0 = z we have

V ij t = e (K-βij I)t v ij + [0,t) e -(K-βij I)(t-s) e p+1 dN j s , = e (K-βij I)t v ij + ξ 1,ij t and Z t = diag(e -ω1t , • • • , e -ω d t )z + t 0 diag(e -ω1(t-s) , • • • , e -ω d (t-s) )dW s , = diag(e -ω1t , • • • , e -ω d t )z + ξ 2,ij t ,
where ξ 1,ij t and ξ 2,ij t are two finite random variables. Since f is bounded it is possible to exchange limits and integrals and one has for any fixed t ≥ 0:

lim (v,z)→∞ T (t)f (v, z) = lim (v,z)→∞ E f (e (K-βij I)t v ij + ξ 1,ij t ) i,j=1,••• ,d , diag(e -ω1t , • • • , e -ω d t )z + ξ 2,ij t , = E lim (v,z)→∞ f (e (K-βij I)t v ij + ξ 1,ij t ) i,j=1,••• ,d , diag(e -ω1t , • • • , e -ω d t )z + ξ 2,ij t , = E [0] , = 0.
This, combined with the weak continuity (trivial) yields that T is a Feller semi-group (cf. the first chapter in [START_REF] Böttcher | Lecture Notes in Mathematics. Lévy-type processes: construction, approximation and sample path properties[END_REF]). Proof. The density of S (E) in Ĉ(E) is a consequence of the density of C ∞ c (E). Let f ∈ S(E) and t ≥ 0. Then

T (t)f (v, z) = E [f (V t , Z t )|(V 0 , Z 0 ) = (v, z)] , = E f (e (K-βij I)t v ij + ξ 1,ij t ) i,j=1,••• ,d , diag(e -ω1t , • • • , e -ω d t )z + ξ 2,ij t , where ξ 1,ij t = [0,t) e -(K-βij I)(t-s) e p+1 dN j s , and ξ 2,ij t = t 0 diag(e -ω1(t-s) , • • • , e -ω d (t-s)
)dW s are two finite random variables. Since S(E) is stable by composition with an affine transform (cf. [START_REF] Galbis | Composition operators on the Schwartz space[END_REF]) and since the differentiation and the expected value can be exchanged (f is bounded), then T (t)f ∈ S(E). Figure 1: The process is self-exciting for α = 2. This is seen as events tend to be clustered in time. Notice how one bin contains two events.

Theorem 3 . 1 .

 31 Let Γ a d-dimension Poisson measure on R 2 + and W a standard d-variate Brownian motion. Let ϕ and δ be functions as defined in the beginning of this subsection. Let µ

2 becomes Assumption 3 . 3 .

 33 The functions ϕ i are l i -Lipschitz for i = 1, • • • , d and the matrix

  We recall that ãm (v, z) = µ m + σ m z m + d j=1 ⟨α mj , v mj ⟩. The jump rates are ϕ 1 (x) = (x) + ∧ 40, ϕ 2 (x) = log(1 + e x ) ∧ 40 and the memory kernels are δ 11 (s) δ 12 (s) δ 21 (s) δ 22 (s) = 5(s -1)e -5s s( s 2 -6)e -6s 0 s( 3s 2 + 4)e -4s .

  ϕ(x) = (x) + ∧ 40 ln(1 + e x ) ∧ 40 4 1+e -(x-3) , and baseline intensities µ = 0.4 1 2 .

  ) = (x) + ∧ 40 ln(1 + e x )

Lemma 7 . 5 .

 75 Let E = (R (p+1) + ) d×d × R d . Let S (E) be the space of Schwartz functions, i.e. the C ∞ functions such that sup x∈E |x m D n f (x)| < +∞, for any multi-indices m and n. Then S (E) is a core for A.

Figure 2 :

 2 Figure2: The process is self-inhibiting for α = -2. Events tend to be distanced from one another.

Figure 3 :

 3 Figure3: A plot of a realisation of a k = (ã m (U k , Y k )) m=1,2 (upper plot) as well as ϕ(a k ) (lower plot) for k ≤ 75. The interaction parameters are (α11 1 , α 11 0 ) = (0, 2), (α 12 1 , α 12 0 ) = (15, -3), (α 21 1 , α 21 0 ) = (-6, 4) and (α22 1 , α 22 0 ) = (10, 0). The perturbation parameters are (σ 1 , σ 2 ) = (0.1, 0.2).

Figure 4 :

 4 Figure 4: A realisation of H k for k ≤ 75.

Figure 5 :

 5 Figure 5: Simulation for h = 0.125. The baseline intensities are (µ 1 , µ 2 ) = (5, 3).

Figure 6 :

 6 Figure 6: Simulation for h = 5 • 10 -3 . The Gaussian perturbation parameters are σ 1 = 1.5 and σ 2 = 0.5. The graph is drawn in simple continuous lines for visibility.

Figure 7 :

 7 Figure 7: The blue curves are the actual kernels and the orange curves are the calibration results. The estimated baseline intensities are μ = (0.37, 1.3, 1.5).

Figure 8 :Figure 9 :

 89 Figure 8: Calibrations of the coefficient α ij 0 and α ij 1 for i, j = 1, • • • , d. The estimation is clearly less precise for α 1 with a larger variance. The red stars are the markers of the ground truth value.

Figure 10 :

 10 Figure 10: The estimator tends to yield larger values for the baseline intensities. This is due to the aforementioned bias, where some self-excitations are ignored.

Figure 11 :Figure 12 :

 1112 Figure 11: The overall effect (excitation, inhibition, independence) are well captured by the exponential kernels. For composite kernels (e.g. δ 12 and δ 22 ) the exponential kernel captures the more influential behaviour.

  Theorem 4.5. Let (U n , Y n , H n ) n∈N be the Markov chain defined in 2.4. ) t∈R+ the càdlàg process that coincides with (U n , Y n , H n ) on the n-th time bin, i.e. ) t∈R+ =⇒ (V t , Z t , N t ) t∈R+ , in distribution as the bin size h goes to zero. Zh t ) t∈R+ to (V t , Z t ) t∈R+ as h -→ 0 (cf. Remark 2.5 and 3.6). Since the Schwartz space S(E) is a core for the infinitesimal generator A (cf. Lemma 7.5), Proposition 4.2 ensures that lim Y ⌊t/h⌋ ), converges to the auxiliary process (V , Z) of the non linear multivariate Hawkes process, weakly in the Skorokhod metric.

	For a bin size h > 0, set (	Ṽ h t ,	Zh t , t ( Ñ h Ṽ h t ,	Zh t ,	Ñ h t ) = (U ⌊t/h⌋ , Y ⌊t/h⌋ , H ⌊t/h⌋ ),
	(cf. Figure 5 for an illustration).
	Then t Proof. It is sufficient to prove the convergence of ( ( Ṽ h t , Zh t , Ñ h	Ṽ h t ,
	Theorem 4.1 yields that the process
				(	Ṽ h t ,	Zh
	To illustrate the convergence, we simulate ãm (	Ṽ h t ,	Zh

h↓0 T ⌊t/h⌋ h f = T (t)f,

for each f in the Banach space (cf. Lemma 7.1) Ĉ(E). t ) = (U ⌊t/h⌋ , t ) in the bi-variate case for t ∈ [0, 15] for two different bin sizes.

Data availability statement

Data sharing not applicable to this article as no datasets were analysed during the current study.

Appendix: Figures


Antoine Voirand
email: a.voirand@brgm.fr

Charles Maragna

Mathilde Veyron

Marc Clausse

A New Borefield Dymola Model Accounting for Groundwater Flow and Complex Borefield Geometries

Keywords: Ground Source Heat Pump, Borehole Heat Exchangers, Borehole Termal Energy Storage, Dymola

In order to support the development of shallow geothermal energy and geothermal energy storage, accurate modelling of ground heat transfer is key. This is particularly true for borehole heat transfer, where the dynamical response of the borehole can greatly affect the overall performance of the surface energy system. Project ANR Sunstone, funded by the French National Agency for Research, aim to optimise the integration of solar heat production and geothermal energy storage, with a practical application on the district heating network of Chateaubriant, Brittany. In this context, the French national survey, BRGM, developed a new borefield Dymola model to contribute to the accuracy of borehole heat transfer modelling. The new borefield model allows considering groundwater flow in the borefield, and heat loss at the surface. It simplifies the declaration of borefield configuration, and allows choosing between rectangular shape on a square base, and hexagonal shape on a hexagonal base. The dynamic heat transfer occurring in a borefield is a complex dynamic process, with time scales involved varying by several order of magnitudes, from the minutes to several decades. The first section of this paper reports the development of new g-functions describing the heat transfer in the borehole, while the second section is focused on its practical implementation in Modelica.

BOREFIELD HEAT TRANSFER MODELLING

Borefield modelling approaches

Heat transfer in the borefield can be split into three main domains of interest. Firstly, the transfer inside the borehole, secondly, the transfer in the ground locally, around each borehole taken individually, and thirdly the thermal interactions between the boreholes inside the borefield at the global level.

Leaving apart the fully discretized models, the transient processes inside the borehole are modelled via an equivalent Thermal Resistance Model (TRM) or Thermal Resistance and Capacity Model (TRCM) [START_REF] Bauer | Thermal Resistance and Capacity Models for Borehole Heat Exchangers[END_REF]), depending on the characteristic time focus of the developed model. TRM is used to simplify the modelling of the short-term behaviour and focus on the long-term processes, while TRCM take into account the influence of the grout thermal capacity and its ability to model the short-term heat transfer, at the expense of the computation time.

The local heat transfer in the ground can be modelled by a discretization of the domain and a corresponding thermal resistance and capacity model. Another approach has been introduced by [START_REF] Eskilson | Thermal Analysis of Heat Extraction Boreholes[END_REF], using a so-called g-function representing the step response of the ground to a constant heat source. The temperature of the ground can be expressed as follows:

𝑇 -𝑇 0 = 𝑞 𝑙 𝜆 𝑔 . 𝑔 [1]
Where T0 is the initial temperature of the ground, ql is the constant linear heat input in the ground, and λg the thermal conductivity of the ground. The step response g can have different expression depending on the geometry of the considered heat source. The simplest one is an Infinite Line Source (ILS) exchanging a constant linear heat flux ql with the ground. However, the ILS model is suitable only for intermediate time.

For short-term behaviour modelling close to the borehole, where the cylindrical nature of the source must be taken into consideration, the Infinite Cylindrical Source (ICS) can be used, and for long-term and long-distance behaviour where the finite nature of the borehole has an impact on the results, the Finite Line Source (FLS) is more suitable. A comprehensive comparison of these basic g-functions and their respective validity domain can be found in [START_REF] Philippe | Validity Ranges of Three Analytical Solutions to Heat Transfer in the Vicinity of Single Boreholes[END_REF].

A wide variety of g-functions has been developed in order to take into account other geometries or boundary conditions, including groundwater flow and boundary condition at the ground surface [START_REF] Conti | Dimensionless Maps for the Validity of Analytical Ground Heat Transfer Models for GSHP Applications[END_REF]). The gfunctions are evaluated before any dynamical modelling, and aggregation methods have been developed in order to use this step-response in the case of varying heat flux input between the borehole and the ground [START_REF] Laferrière | Development and Validation of a Full-Time-Scale Semi-Analytical Model for the Short-and Long-Term Simulation of Vertical Geothermal Bore Fields[END_REF].

Some g-functions have an analytic expression, but as the modelling hypothesis become more and more general, the difficulty of finding such an expression becomes important.

At the global level, g-functions can also be used to compute the thermal interactions between boreholes. Another approach, generally implemented when the Borefield is used as a Borehole Thermal Energy Storage (BTES), consists of a discretization of the ground inside the storage, and the computation of an average ground temperature for each sub-domain, based on a local modelling of a borehole heat transfer and the interaction with the adjacent sub-domain. This approach proceeds by a wide approximation of the heat transfer processes in the ground but can be very timesaving and sufficiently precise if carefully parameterized. Furthermore, it allows the calculation of the mean temperature of the borefield ground, which is a valuable information for the control of the charge and discharge cycles of the borefield when taken as a heat storage medium.

Borefield modelling in Modelica language

The Modelica language is particularly efficient for complex equation based multi-physics modelling and simulation. It supports acausal connection and objectoriented construction of components and is therefore adapted for complex energy systems simulation such as the ground source heat pump dynamic modelling and borefield thermal performances optimisation.

Two borefield models have already been implemented in the Modelica language. The first and more widely used is the Hybrid Step Response Model (HSRM) developed by [START_REF] Picard | Advanced Hybrid Model for Borefield Heat Exchanger Performance Evaluation, an Implementation in Modelica[END_REF], and included in the IBPSA and BUILDINGS model libraries. It uses a TRCM model for short-time accuracy, and an asymptotic matching between ICS and FLS g-functions for the local heat transfer in the ground [START_REF] Li | Full-Scale Temperature Response Function (G-Function) for Heat Transfer by Borehole Ground Heat Exchangers (GHEs) from Sub-Hour to Decades[END_REF], allowing the resulting g-function to be accurate at both short-term and long-term delay. Interactions between boreholes are calculated using only the FLS method, since the boreholes are separated by several meters. An aggregation method is used to evaluate varying heat flux in a timesaving way, but the calculation of the gfunction step response during the initialisation can be time consuming for large borefield.

Another model is the MOBTES library developed by [START_REF] Formhals | A Modelica Toolbox for the Simulation of Borehole Thermal Energy Storage Systems[END_REF]. Its main difference is the global level which is based on the discretisation approach in the same way that the Duct Ground Storage model (DST) implemented in the dynamic modelling software TRNSYS and developed by Pahud and Hellström (1996).

These two models do not take into account the influence of groundwater flow. [START_REF] Formhals | A Modelica Toolbox for the Simulation of Borehole Thermal Energy Storage Systems[END_REF] justify this omission by the fact that BTES are generally installed in areas where there is no groundwater flow, since its impact is largely negative on the BTES thermal performances. That is indeed true, but the MOBTES model is therefore not capable of modelling the actual impact of groundwater flow and considering possible mitigation actions. There are also cases where the borefield is not considered as a thermal energy storage system but rather as a heating and cooling system connected to Heat Pump, that can be installed in areas where groundwater flow is present. It is therefore important to be able to consider the effect on the borefield thermal performance.

It is possible to analytically derive a model of Moving Finite Line Source (MFLS) or Moving Infinite Line Source (MILS) from the FLS and ILS models respectively. It is however impossible to do so with a Moving Infinite Cylindrical Source (MICS) from the ICS model. [START_REF] Conti | Transient Forced Convection from an Infinite Cylindrical Heat Source in a Saturated Darcian Porous Medium[END_REF] perform a numerical study based on the same assumptions than the MICS and give g-functions approximation for the temperature difference from the initial temperature at the borehole surface for discrete values of Peclet number. However, discrete g-functions do not allow for accurate representation of the temperature step response over the entire spectrum of groundwater velocities, and the MICS model is accurate for shortterm modelling only since the influence of the finite length of the borehole can have major impact on the step response for long-term modelling.

Numerical calculation of a Moving Finite Cylindrical Heat Source model

The Finite Element based simulation software COMSOL Multiphysics v4.2a was used to evaluate the temperature step response at the borehole wall of a finite cylindrical heat source with groundwater flow and mixed boundary condition at the surface (MFCSm). The characteristic length of the problem is the borehole radius rb. The Fourier number, the nondimensional borehole depth, the Peclet number and the Biot number between the convective heat transfer coefficient h at the surface and the ground conductivity λg are defined as follows:

𝐹𝑜 = 𝛼 𝑒𝑞 𝑡 𝑟 𝑏 [2] 𝐻 * = 𝐻 𝑟 𝑏 [3] 𝑃𝑒 = 𝑢 𝑇 𝑟 𝑏 (𝜌𝑐 𝑝 ) 𝑔 𝜆 𝑔 [4] 𝐵𝑖 = ℎ 𝑟 𝑏 𝜆 𝑔 [5]
With αeq the equivalent heat diffusivity of the waterfilled ground, function of the porosity of the ground, H the borehole length and uT the effective heat transport velocity, defined from the Darcy velocity u in the principal direction of the flow:

𝑢 𝑇 = 𝑢 (𝜌𝑐 𝑝 ) 𝑤 (𝜌𝑐 𝑝 ) 𝑔 [6]
Where the subscript w refers to the physical quantities for water and the subscript g refers to the ground.

Figure 1 shows the influence of the Peclet number on the g-function step-response at the borehole wall for the numerical solution from COMSOL, the analytical solution MFLS with imposed temperature at the surface (i.e. 𝐵𝑖 → ∞), and the ICS solution with no groundwater flow. Groundwater flow has no impact on the step-response at low Fourier number but leads to the establishment of asymptotes for high Fo. Figure 1 shows a good agreement between ICS and the COMSOL solution for low Fo, and good agreement between COMSOL and MFLS for high Fo number. The COMSOL model is thus a good approximation for both short time and long-time calculation of the stepresponse of the borehole wall.

Figure 2 displays the g-function step-response of the borehole wall temperature for different values of Biot number with no groundwater flow and H * = 100, from the COMSOL numerical solution and from a MFLS model with mixed thermal conditions at surface (MFLSm). The analytical, quite complex, expressions of this later model can be found in [START_REF] Rivera | A Finite Line Source Model with Cauchy-Type Top Boundary Conditions for Simulating near Surface Effects on Borehole Heat Exchangers[END_REF]. Figure 2 shows a good agreement between the numerical and analytical solutions for high Fourier numbers. However, the numerical solution is more accurate at low Fourier number. The meta-model is validated by comparing the COMSOL solution, the fitted curve on this solution, and the result of the meta-model for several parameter sets not used for the elaboration of the meta-model. Figure 3 shows an example of this validation.

Meta-model derivation of the MFCSm model from the COMSOL numerical model

The draft of a review article presenting details about the COMSOL model, and the elaboration of the metamodel is currently under way.

IMPLEMENTATION OF THE MOVING FINITE CYLINDRICAL SOURCE

The meta-model derived in section 1.4 is implemented in a new Modelica library. This library is based on the same structure than the HSRM model developed by [START_REF] Picard | Advanced Hybrid Model for Borefield Heat Exchanger Performance Evaluation, an Implementation in Modelica[END_REF], since the calculated function is the borehole wall temperature response to a step of power in both cases. The MICSm meta-model calculates the g-function for a single borehole in the borefield, and the MFLS model calculates the interactions between boreholes, since the cylindrical nature of the heat source is negligible at such distances. This modification allows for the consideration of groundwater flow in the borefield, and mixed boundary conditions at the surface.

At the difference with the HSRM model, only one model is proposed at the top level of the arborescence (Figure 4). The user can then choose in the parameter tab of this model (Figure 5) the type of borehole used between three choices: single-U, parallel double-U or series double-U. It is also possible to choose if there is a groundwater flow. If no groundwater flow is chosen, the HSRM model is used. Moreover, the borefield configuration definition in the HSRM model can only be achieved by providing the cartesian coordinates of all the boreholes, which can be time consuming and cumbersome, or require elaborate function declaration in the configuration record. The new borefield Modelica model presented here does not change the calculation process of borehole interactions, so the calculus core still requires the coordinates of the entire borefield, but it simplifies the declaration process and give access to two different patterns on which the boreholes are placed: rectangular borefield on a square base, or hexagonal borefield on a hexagonal base (Figure 6). The first configuration requires the following parameters: number of rows (x-direction), number of columns (y-direction) in the borefield, and distance between the boreholes. By default, the direction of groundwater flow is parallel to the rows. It is to this date not possible to specify a direction of groundwater flow with respect to the row direction. The second configuration requires the following parameters: number of boreholes, and distance between boreholes. The borefield is then constructed by progressively filling the concentric annular, on a hexagonal base, around a central borehole.

CONCLUSIONS

The Modelica language is particularly adapted for dynamic multi-physics energy system simulation. Shallow geothermal energy systems need to rely on accurate modelling of the heat transfer occurring between the borehole and the ground. However, the available borehole and borefield Modelica models in the commonly used libraries do not consider the impact of groundwater flow on the temperature step-response of the borehole or the borefield.

The new borefield Modelica module presented here considers groundwater flow and mixed boundary conditions. It is based on a meta-model derived from COMSOL simulations of the temperature step-response at the borehole wall function of the non-dimensional borehole height, the Peclet number and the Biot numbers.

Four

  hundred and thirty numerical evaluations of the gfunction step-response of the MFCSm model with parameters chosen using a Latin hypercube sample in the following ranges: step-responses were approximated using modified sigmoids: 𝐺 𝑀𝐹𝐶𝑆 𝑚 (𝐻 * , 𝑃𝑒, 𝐵𝑖, 𝐹𝑜) (𝐹𝑜)) 𝑛 )) [10] 𝑌 = [log 10 𝐻 * log 10 𝐵𝑖 log 10 𝑃𝑒] [11] The coefficients Bn (n = 0, 1, 2, 3) and the asymptote Gmax are computed as a polynomial function of Y: vectors 𝑎 𝐵 𝑛 and 𝑎 𝐺 𝑚𝑎𝑥 contain 35 coefficients. [𝑚 𝑖𝑗 ] has 35 lines and 3 columns: for instance, a line with [0 1 2] corresponds to ∏ 𝑌 𝑖 𝑚 equation [12] and [13]. The coefficients are found numerically by minimising the standard deviation between the g-function computed by COMSOL and eq. 10.

Figure

  Figure 1: ICS model (black) vs COMSOL solution (solid-coloured lines) vs MFLS (dashed) function of the Fourier number for large Biot number (imposed surface temperature) and H * = 100.

Figure

  Figure 2: COMSOL solution (solid) vs MFLSm (dashed) function of the Fourier number for Pe = 0 and H * = 100.

Figure 3 :

 3 Figure 3: COMSOL solution of the temperature wall step-response vs fitted function based on eq (10) vs correlated function from the metamodel. H * = 320; Pe = 0.03; Bi = 0.48.

Figure 4 :

 4 Figure 4: Arborescence of the new Modelica borefield module.

Figure 5 :

 5 Figure 5: Parameter display of the borefield model.

Figure 6 :

 6 Figure 6: Parameter display of the configuration data record in the new Modelica borefield module.

Acknowledgements

The developments presented in this article have been funded by the French National Research Agency under the grant ANR-17-CE05-0035 (SUNSTONE)