Benoît Sagot
email: benoit.sagot@inria.fr

Éric De La Clergerie
email: eric.de_la_clergerie@inria.fr

Error mining in parsing results

We introduce an error mining technique for automatically detecting errors in resources that are used in parsing systems. We applied this technique on parsing results produced on several million words by two distinct parsing systems, which share the syntactic lexicon and the pre-parsing processing chain. We were thus able to identify missing and erroneous information in these resources.

Introduction

Natural language parsing is a hard task, partly because of the complexity and the volume of information that have to be taken into account about words and syntactic constructions. However, it is necessary to have access to such information, stored in resources such as lexica and grammars, and to try and minimize the amount of missing and erroneous information in these resources. To achieve this, the use of these resources at a largescale in parsers is a very promising approach [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF], and in particular the analysis of situations that lead to a parsing failure: one can learn from one's own mistakes.

We introduce a probabilistic model that allows to identify forms and form bigrams that may be the source of errors, thanks to a corpus of parsed sentences. In order to facilitate the exploitation of forms and form bigrams detected by the model, and in particular to identify causes of errors, we have developed a visualization environment. The whole system has been tested on parsing results produced for several multi-million-word corpora and with two different parsers for French, namely SXLFG and FRMG.

However, the error mining technique which is the topic of this paper is fully system-and language-independent. It could be applied without any change on parsing results produced by any system working on any language. The only information that is needed is a boolean value for each sentence which indicates if it has been successfully parsed or not.

Principles

General idea

The idea we implemented is inspired from [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF]. In order to identify missing and erroneous information in a parsing system, one can analyze a large corpus and study with statistical tools what differentiates sentences for which parsing succeeded from sentences for which it failed.

The simplest application of this idea is to look for forms, called suspicious forms, that are found more frequently in sentences that could not be parsed. This is what van [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF] does, without trying to identify a suspicious form in any sentence whose parsing failed, and thus without taking into account the fact that there is (at least) one cause of error in each unparsable sentence. 1On the contrary, we will look, in each sentence on which parsing failed, for the form that has the highest probability of being the cause of this failure: it is the main suspect of the sentence. This form may be incorrectly or only partially described in the lexicon, it may take part in constructions that are not described in the grammar, or it may exemplify imperfections of the pre-syntactic processing chain. This idea can be easily extended to sequences of forms, which is what we do by tak-ing form bigrams into account, but also to lemmas (or sequences of lemmas).

Form-level probabilistic model

We suppose that the corpus is split in sentences, sentences being segmented in forms. We denote by s i the i-th sentence. We denote by o i,j , (1 ≤ j ≤ |s i |) the occurrences of forms that constitute s i , and by F (o i,j ) the corresponding forms. Finally, we call error the function that associates to each sentence s i either 1, if s i 's parsing failed, and 0 if it succeeded.

Let O f be the set of the occurrences of a form f in the corpus:

O f = {o i,j |F (o i,j ) = f }. The number of occurrences of f in the corpus is there- fore |O f |.
Let us define at first the mean global suspicion rate S, that is the mean probability that a given occurrence of a form be the cause of a parsing failure. We make the assumption that the failure of the parsing of a sentence has a unique cause (here, a unique form. . . ). This assumption, which is not necessarily exactly verified, simplifies the model and leads to good results. If we call occ total the total amount of forms in the corpus, we have then:

S = Σ i error(s i ) occ total
Let f be a form, that occurs as the j-th form of sentence s i , which means that F (o i,j ) = f . Let us assume that s i 's parsing failed: error(s i ) = 1. We call suspicion rate of the j-th form o i,j of sentence s i the probability, denoted by S i,j , that the occurrence o i,j of form form f be the cause of the s i 's parsing failure. If, on the contrary, s i 's parsing succeeded, its occurrences have a suspicion rate that is equal to zero.

We then define the mean suspicion rate S f of a form f as the mean of all suspicion rates of its occurrences:

S f = 1 |O f | • o i,j ∈O f S i,j
To compute these rates, we use a fix-point algorithm by iterating a certain amount of times the following computations. Let us assume that we just completed the n-th iteration: we know, for each sentence s i , and for each occurrence o i,j of this sentence, the estimation of its suspicion rate S i,j as computed by the n-th iteration, estimation that is denoted by S (n) i,j . From this estimation, we compute the n + 1-th estimation of the mean suspicion rate of each form f , denoted by S (n+1) f :

S (n+1) f = 1 |O f | • o i,j ∈O f S (n) i,j
This rate 2 allows us to compute a new estimation of the suspicion rate of all occurrences, by giving to each occurrence if a sentence s i a suspicion rate S (n+1) i,j that is exactly the estimation S (n+1) f of the mean suspicion rate of S f of the corresponding form, and then to perform a sentencelevel normalization. Thus:

S (n+1) i,j = error(s i ) • S (n+1) F (o i,j ) 1≤j≤|s i | S (n+1) F (o i,j )
At this point, the n+1-th iteration is completed, and we can resume again these computations, until convergence on a fix-point. To begin the whole process, we just say, for an occurrence o i,j of sentence s i , that S (0) i,j = error(s i )/|s i |. This means that for a non-parsable sentence, we start from a baseline where all of its occurrences have an equal probability of being the cause of the failure.

After a few dozens of iterations, we get stabilized estimations of the mean suspicion rate each form, which allows:

• to identify the forms that most probably cause errors,

• for each form f , to identify non-parsable sentences s i where an occurrence o i,j ∈ O f of f is a main suspect and where o i,j has a very 2 We also performed experiment in which S f was estimated by an other estimator, namely the smoothed mean suspicion rate, denoted by S(n) f , that takes into account the number of occurrences of f . Indeed, the confidence we can have in the estimation S (n) f is lower if the number of occurrences of f is lower. Hence the idea to smooth S

-λ depend on |O f |: if |O f | is high, S(n) f will be close from S (n) f
; if it is low, it will be closer from S:

S(n) f = λ(|O f |) • S (n) f + (1 -λ(|O f |)) • S.
In these experiments, we used the smoothing function

λ(|O f |) = 1 -e -β|O f | with β = 0.1. But this model, used with the ranking according to M f = S f • ln |O f | (see below)
, leads results that are very similar to those obtained without smoothing. Therefore, we describe the smoothingless model, which has the advantage not to use an empirically chosen smoothing function.

high suspicion rate among all occurrences of form f .

We implemented this algorithm as a perl script, with strong optimizations of data structures so as to reduce memory and time usage. In particular, form-level structures are shared between sentences.

Extensions of the model

This model gives already very good results, as we shall see in section 4. However, it can be extended in different ways, some of which we already implemented.

First of all, it is possible not to stick to forms. Indeed, we do not only work on forms, but on couples made out of a form (a lexical entry) and one or several token(s) that correspond to this form in the raw text (a token is a portion of text delimited by spaces or punctuation tokens).

Moreover, one can look for the cause of the failure of the parsing of a sentence not only in the presence of a form in this sentence, but also in the presence of a bigram3 of forms. To perform this, one just needs to extend the notions of form and occurrence, by saying that a (generalized) form is a unigram or a bigram of forms, and that a (generalized) occurrence is an occurrence of a generalized form, i.e., an occurrence of a unigram or a bigram of forms. The results we present in section 4 includes this extension, as well as the previous one.

Another possible generalization would be to take into account facts about the sentence that are not simultaneous (such as form unigrams and form bigrams) but mutually exclusive, and that must therefore be probabilized as well. We have not yet implemented such a mechanism, but it would be very interesting, because it would allow to go beyond forms or n-grams of forms, and to manipulate also lemmas (since a given form has usually several possible lemmas).

Experiments

In order to validate our approach, we applied these principles to look for error causes in parsing results given by two deep parsing systems for French, FRMG and SXLFG, on large corpora.

Parsers

Both parsing systems we used are based on deep non-probabilistic parsers. They share:

• the Lefff 2 syntactic lexicon for French (Sagot et al., 2005), that contains 500,000 entries (representing 400,000 different forms) ; each lexical entry contains morphological information, sub-categorization frames (when relevant), and complementary syntactic information, in particular for verbal forms (controls, attributives, impersonals,. . . ),

• the SXPipe pre-syntactic processing chain (Sagot and Boullier, 2005), that converts a raw text in a sequence of DAGs of forms that are present in the Lefff ; SXPipe contains, among other modules, a sentence-level segmenter, a tokenization and spelling-error correction module, named-entities recognizers, and a non-deterministic multi-word identifier.

But FRMG and SXLFG use completely different parsers, that rely on different formalisms, on different grammars and on different parser builder. Therefore, the comparison of error mining results on the output of these two systems makes it possible to distinguish errors coming from the Lefff or from SXPipe from those coming to one grammar or the other. Let us describe in more details the characteristics of these two parsers.

The FRMG parser (Thomasset and Villemonte de la Clergerie, 2005) is based on a compact TAG for French that is automatically generated from a meta-grammar. The compilation and execution of the parser is performed in the framework of the DYALOG system (Villemonte de la [START_REF] Villemonte De La Clergerie | DyALog: a tabular logic programming based environment for NLP[END_REF].

The SXLFG parser (Boullier and Sagot, 2005b;Boullier and Sagot, 2005a) is an efficient and robust LFG parser. Parsing is performed in two steps. First, an Earley-like parser builds a shared forest that represents all constituent structures that satisfy the context-free skeleton of the grammar. Then functional structures are built, in one or more bottom-up passes. Parsing efficiency is achieved thanks to several techniques such as compact data representation, systematic use of structure and computation sharing, lazy evaluation and heuristic and almost non-destructive pruning during parsing.

Both parsers implement also advanced error recovery and tolerance techniques, but they were useless for the experiments described here, since we want only to distinguish sentences that receive a full parse (without any recovery technique) from those that do not.

Corpora

We parsed with these two systems the following corpora:

MD corpus : This corpus is made out of 14.5 million words (570,000 sentences) of general journalistic corpus that are articles from the Monde diplomatique.

EASy corpus : This is the 40,000-sentence corpus that has been built for the EASy parsing evaluation campaign for French [START_REF] Paroubek | EASy : campagne d'évaluation des analyseurs syntaxiques[END_REF]. We only used the raw corpus (without taking into account the fact that a manual parse is available for 10% of all sentences). The EASy corpus contains several sub-corpora of varied style: journalistic, literacy, legal, medical, transcription of oral, email, questions, etc.

Both corpora are raw in the sense that no cleaning whatsoever has been performed so as to eliminate some sequences of characters that can not really be considered as sentences.

Table 1 gives some general information on these corpora as well as the results we got with both parsing systems. It shall be noticed that both parsers did not parse exactly the same set and the same number of sentences for the MD corpus, and that they do not define in the exactly same way the notion of sentence.

Results visualization environment

We developed a visualization tool for the results of the error mining, that allows to examine and annotate them. It has the form of an HTML page that uses dynamic generation methods, in particular javascript. An example is shown on Figure 1.

To achieve this, suspicious forms are ranked according to a measure M f that models, for a given form f , the benefit there is to try and correct the (potential) corresponding error in the resources. A user who wants to concentrate on almost certain errors rather than on most frequent ones can visualize suspicious forms ranked according to M f = S f . On the contrary, a user who wants to concentrate on most frequent potential errors, rather than on the confidence that the algorithm has given to errors, can visualize suspicious forms ranked according to 4 M f = S f |O f |. The default choice, which is adopted to produce all tables shown in this paper, is a balance between these two possibilities, and ranks suspicious forms according to

M f = S f • ln |O f |.
The visualization environment allows to browse through (ranked) suspicious forms in a scrolling list on the left part of the page (A). When the suspicious form is associated to a token that is the same as the form, only the form is shown. Otherwise, the token is separated from the form by the symbol " / ". The right part of the page shows various pieces of information about the currently selected form. After having given its rank according to the ranking measure M f that has been chosen (B), a field is available to add or edit an annotation associated with the suspicious form (D). These annotations, aimed to ease the analysis of the error mining results by linguists and by the developers of parsers and resources (lexica, grammars), are saved in a database (SQLITE). Statistical information is also given about f (E), including its number of occurrences occ f , the number of occurrences of f in non-parsable sentences, the final estimation of its mean suspicion rate S f and the rate err(f ) of non-parsable sentences among those where f appears. This indications are complemented by a brief summary of the iterative process that shows the convergence of the successive estimations of S f . The lower part of the page gives a mean to identify the cause of f -related errors by showing 4 Let f be a form. The suspicion rate S f can be considered as the probability for a particular occurrence of f to cause a parsing error. Therefore, S f |O f | models the number of occurrences of f that do cause a parsing error. f 's entries in the Lefff lexicon (G) as well as nonparsable sentences where f is the main suspect and where one of its occurrences has a particularly high suspicion rate 5 (H).

The whole page (with annotations) can be sent by e-mail, for example to the developer of the lexicon or to the developer of one parser or the other (C).

Results

In this section, we mostly focus on the results of our error mining algorithm on the parsing results provided by SXLFG on the MD corpus. We first present results when only forms are taken into account, and then give an insight on results when both forms and form bigrams are considered.

5 Such an information, which is extremely valuable for the developers of the resources, can not be obtained by global (form-level and not occurrence-level) approaches such as the err(f )-based approach of [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF]. Indeed, enumerating all sentences which include a given form f , and which did not receive a full parse, is not precise enough: it would show at the same time sentences wich fail because of f (e.g., because its lexical entry lacks a given subcategorization frame) and sentences which fail for an other independent reason.

Finding suspicious forms

The execution of our error mining script on MD/SXLFG, with i max = 50 iterations and when only (isolated) forms are taken into account, takes less than one hour on a 3.2 GHz PC running Linux with a 1.5 Go RAM. It outputs 18,334 relevant suspicious forms (out of the 327,785 possible ones), where a relevant suspicious form is defined as a form f that satisfies the following arbitrary constraints:6 S (imax) f > 1, 5 • S and |O f | > 5. We still can not prove theoretically the convergence of the algorithm. 7 But among the 1000 bestranked forms, the last iteration induces a mean variation of the suspicion rate that is less than 0.01%.

On a smaller corpus like the EASy corpus, 200 iterations take 260s. The algorithm outputs less than 3,000 relevant suspicious forms (out of the 61,125 possible ones). Convergence information is the same as what has been said above for the MD corpus.

Table 2 gives an idea of the repartition of suspicious forms w.r.t. their frequency (for FRMG on MD), showing that rare forms have a greater probability to be suspicious. The most frequent suspicious form is the double-quote, with (only) S f = 9%, partly because of segmentation problems.

Analyzing results

Table 3 gives an insight on the output of our algorithm on parsing results obtained by SXLFG on the MD corpus. For each form f (in fact, for each couple of the form (token,form)), this table displays its suspicion rate and its number of occurrences, as well as the rate err(f ) of non-parsable sentences among those where f appears and a short manual analysis of the underlying error.

In fact, a more in-depth manual analysis of the results shows that they are very good: errors are correctly identified, that can be associated with four error sources: (1) the Lefff lexicon, (2) the SXPipe pre-syntactic processing chain, (3) imperfections of the grammar, but also (4) problems related to the corpus itself (and to the fact that it is a raw corpus, with meta-data and typographic noise).

On the EASy corpus, results are also relevant, but sometimes more difficult to interpret, because of the relative small size of the corpus and because of its heterogeneity. In particular, it contains email and oral transcriptions sub-corpora that introduce a lot of noise. Segmentation problems (caused both by SXPipe and by the corpus itself, which is already segmented) play an especially important role.

Comparing results with results of other algorithms

In order to validate our approach, we compared our results with results given by two other relevant algorithms:

• van Noord's [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF]) (form-level and non-iterative) evaluation of err(f ) (the rate of non-parsable sentences among sentences containing the form f ),

• a standard (occurrence-level and iterative) maximum entropy evaluation of each form's contribution to the success or the failure of a sentence (we used the MEGAM package [START_REF] Daumé | Notes on CG and LM-BFGS optimization of logistic regression[END_REF]).

As done for our algorithm, we do not rank forms directly according to the suspicion rate S f computed by these algorithms. Instead, we use the M f measure presented above (M f = S f •ln |O f |). Using directly van Noord's measure selects as most suspicious words very rare words, which shows the importance of a good balance between suspicion rate and frequency (as noted by [START_REF] Van Noord | Error mining for widecoverage grammar engineering[END_REF] in the discussion of his results). This remark applies to the maximum entropy measure as well.

Table 4 shows for all algorithms the 10 bestranked suspicious forms, complemented by a manual evaluation of their relevance. One clearly sees that our approach leads to the best results. Van Noord's technique has been initially designed to find errors in resources that already ensured a very high coverage. On our systems, whose development is less advanced, this technique ranks as most suspicious forms those which are simply the most frequent ones. It seems to be the case for the standard maximum entropy algorithm, thus showing the importance to take into account the fact that there is at least one cause of error in any sentence whose parsing failed, not only to identify a main suspicious form in each sentence, but also to get relevant global results.

Comparing results for both parsers

We complemented the separated study of error mining results on the output of both parsers by an analysis of merged results. We computed for each form the harmonic mean of both measures M f = S f • ln |O f | obtained for each parsing system. Results (not shown here) are very interesting, because they identify errors that come mostly from resources that are shared by both systems (the Lefff lexicon and the pre-syntactic processing chain SXPipe). Although some errors come from common lacks of coverage in both grammars, it is nevertheless a very efficient mean to get a first repartition between error sources.

Introducing form bigrams

As said before, we also performed experiments where not only forms but also form bigrams are treated as potential causes of errors. This approach allows to identify situations where a form is not in itself a relevant cause of error, but leads often to a parse failure when immediately followed or preceded by an other form. phasize bigram results), with the same data as in table 3.

Conclusions and perspectives

As we have shown, parsing large corpora allows to set up error mining techniques, so as to identify missing and erroneous information in the different resources that are used by full-featured parsing systems. The technique described in this paper and its implementation on forms and form bigrams has already allowed us to detect many errors and omissions in the Lefff lexicon, to point out inappropriate behaviors of the SXPipe pre-syntactic processing chain, and to reveal the lack of coverage of the grammars for certain phenomena.

We intend to carry on and extend this work. First of all, the visualization environment can be enhanced, as is the case for the implementation of the algorithm itself.

We would also like to integrate to the model the possibility that facts taken into account (today, forms and form bigrams) are not necessarily certain, because some of them could be the consequence of an ambiguity. For example, for a given form, several lemmas are often possible. The probabilization of these lemmas would thus allow to look for most suspicious lemmas.

We are already working on a module that will allow not only to detect errors, for example in the lexicon, but also to propose a correction. To achieve this, we want to parse anew all nonparsable sentences, after having replaced their main suspects by a special form that receives under-specified lexical information. These information can be either very general, or can be computed by appropriate generalization patterns applied on the information associated by the lexicon with the original form. A statistical study of the new parsing results will make it possible to propose corrections concerning the involved forms.

Figure 1 :

 1 Figure 1: Error mining results visualization environment (results are shown for MD/FRMG).

Table 1 :

 1 General information on corpora and parsing results

	corpus	#sentences	#success (%)	#forms	#occ	S (%)	Date
	MD/FRMG	330,938	136,885 (41.30%) 255,616 10,422,926 1.86% Jul. 05
	MD/SXLFG	567,039	343,988 (60.66%) 327,785 14,482,059 1.54% Mar. 05
	EASy/FRMG	39,872	16,477 (41.32%)	61,135	878,156	2.66% Dec. 05
	EASy/SXLFG	39,872	21,067 (52.84%)	61,135	878,156	2.15% Dec. 05

Table 5

 5

	shows best-ranked form bigrams (forms
	that are ranked in-between are not shown, to em-

Table 2 :

 2 Suspicious forms repartition for MD/FRMG

	Rank	Token(s)/form	S	(50) f	|O f | err(f ) M f	Error cause
	1	_____/_UNDERSCORE 100% 4804	100%	8.48 corpus: typographic noise
	2	poussé	76%	91	94%	3.41 Lefff : misses as an adjective
	3	(...)	46%	1216	74%	3.27 SXPipe: should be treated as skippable words
	4	fallu	55%	215	93%	2.94 Lefff : should allow nominal object (not only infinitive)
	5	organisés	62%	89	94%	2.78 Lefff : misses as an adjective
	7	puisse	48%	287	95%	2.74 grammar: construction missing (aussi adj que puisse inf. . . )
	8	fallait	51%	221	91%	2.73 as for fallu
	9	décentralisée	55%	128	88%	2.68 Lefff : misses as an adjective
	10	généralisé	56%	115	96%	2.64 Lefff : misses as an adjective
	11	spécialisées	58%	90	94%	2.63 Lefff : misses as an adjective

Table 3 :

 3 Analysis of the 10 best-ranked forms (ranked according toM f = S f • ln |O f |)

		this paper		global		maxent
	Rank	Token(s)/form	Eval	Token(s)/form Eval	Token(s)/form Eval
	1	_____/_UNDERSCORE	++	*	+	pour	-
	2	(...)	++	,	-	)	-
	3	2_]/_NUMBER	++	livre	-	à	-
	4	privées	++	.	-	qu'il/qu'	-
	5	Haaretz/_Uw	++	de	-	sont	-
	6	contesté	++	;	-	le	-
	7	occupés	++	:	-	qu'un/qu'	+
	8	privée	++	la	-	qu'un/un	+
	9	[...]	++	étrangères	-	que	-
	10	faudrait	++	lecteurs	-	pourrait	-

Table 4 :

 4 The 10 best-ranked suspicious forms, according the the M f measure, as computed by different algorithms: ours (this paper), a standard maximum entropy algorithm (maxent) and van Noord's rate err(f ) (global).

	Rank Tokens and forms M f	Error cause
	4	Toutes/toutes les	2.73 grammar: badly treated pre-determiner adjective
	6	y en	2,34 grammar: problem with the construction il y en a. . .
	7	in "	1.81 Lefff : in misses as a preposition, which happends before book titles (hence the ")
	10	donne à	1.44 Lefff : donner should sub-categorize à-vcomps (donner à voir. . . )
	11	de demain	1.19 Lefff : demain misses as common noun (standard adv are not preceded by prep)
	16	( 22/_NUMBER	0.86 grammar: footnote references not treated
	16	22/_NUMBER )	0.86 as above

Table 5 :

 5 Best ranked form bigrams (forms ranked inbetween are not shown; ranked according to M f = S

f • ln |O f |). These results have been computed on a subset of the MD corpus (60,000 sentences).

Indeed, he defines the suspicion rate of a form f as the rate of unparsable sentences among sentences that contain f .

One could generalize this to n-grams, but as n gets higher the number of occurrences of n-grams gets lower, hence leading to non-significant statistics.

These constraints filter results, but all forms are taken into account during all iterations of the algorithm.

However, the algorithms shares many common points with iterative algorithm that are known to converge and that have been proposed to find maximum entropy probability distributions under a set of constraints[START_REF] Berger | A maximun entropy approach to natural language processing[END_REF]. Such an algorithm is compared to ours later on in this paper.


Alexandre Brouste

Christophe Dutang

Tom Rohmer

Closed-form Maximum Likelihood Estimator for Generalized Linear Models in the case of categorical explanatory variables: Application to insurance loss modeling

Keywords: Regression models, heavy-tailed distributions, explicit MLE, insurance claim modeling

come

Introduction

The assumption of identical distributions for random variables in an observation sample is relaxed for regression models by considering explanatory variables. Generalized Linear Models (GLMs) were introduced by [START_REF] Nelder | Generalized linear models[END_REF] and popularized in [START_REF] Mccullagh | Generalized linear models[END_REF]. GLMs rely on probability distribution functions of exponential type for the response variable which include most of the light and medium tailed distributions (such as normal, gamma or inverse Gaussian). Asymptotic properties of sequences of maximum likelihood estimators (MLE) for GLMs were studied by [START_REF] Fahrmeir | Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models[END_REF].

The finite sample property of MLE of specific GLMs has been discussed in depth in statistical literature, see [START_REF] Fienberg | The Analysis of Cross-Classified Categorical Data[END_REF] and [START_REF] Haberman | Log-linear models for frequency tables with ordered classifications[END_REF]. In addition, there are lots of literature studying the finite sample property of MLE for logistic regression models, see e.g. [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF] and [START_REF] Silvapulle | On the existence of maximum likelihood estimators for the binomial response models[END_REF].

Regression models for heavy-tailed distributions have been mainly studied through the point-of-view of extreme value analysis, see [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF] for a review. A regression model for the generalized Pareto distribution (GPD) where the scale parameter depends on covariates are described in [START_REF] Davison | Models for exceedances over high thresholds[END_REF] with a least square estimation procedure and a model checking method. [START_REF] Beirlant | Burr regression and portfolio segmentation[END_REF] propose a Burr model by regressing the shape parameter with an exponential link on explanatory variables. In the aforementioned article, a simulation study with one explanatory variable is detailed as well as an application to fire insurance. Residual plots and asymptotic convergence towards the normal distribution are also discussed. Similarly, [START_REF] Ozkok | Bayesian modelling of the time delay between diagnosis and settlement for critical illness insurance using a burr generalised-linear-type model[END_REF] propose a regression model for Burr distribution where the scale parameter depends on covariates.

An estimation of the extremal tail index (used in generalized extreme value (GEV) distributions and GPD) by considering a class of distribution function with an exponential link on explanatory variables is also described in [START_REF] Beirlant | Regression with response distributions of pareto-type[END_REF]. Using generalized residuals of explanatory variable makes possible the estimation of the tail index. Still by the extreme value theory approach, [START_REF] Chavez-Demoulin | An extreme value approach for modeling operational risk losses depending on covariates[END_REF] and [START_REF] Hambuckers | A semiparametric model for generalized pareto regression based on a dimension reduction assumption[END_REF] both propose a semi-parametric regression model for GEV and GPD where the explanatory variables are time or known factor levels. They assume that all parameters depend on covariates and also use exponentially distributed residuals.

Outside the extreme value theory framework, there is also a literature studying the covariate modelling. For instance, the so-called double GLMs, where the dispersion parameter also depends on covariates, have been studied by McCullagh and Nelder (1989, Chap. 10) or [START_REF] Smyth | Adjusted likelihood methods for modelling dispersion in generalized linear models[END_REF]. Furthermore, [START_REF] Rigby | Generalized additive models for location, scale and shape[END_REF] propose a general regression framework where all parameters are modeled by explanatory variable and the distribution is not restricted to exponential family. The only restriction that the authors impose is the twice differentiability of the density function w.r.t. parameters. However, there is no clear convergence result of the proposed estimators. Among the proposed distributions, [START_REF] Rigby | Generalized additive models for location, scale and shape[END_REF] use 1-parameter Pareto, log-logistic (a special case of the Pareto 3 distribution) and GEV distributions.

In this paper, we propose closed-form estimators for GLMs in the case of categorical variables. The expression is valid for any distribution belonging to the one-parameter exponential family and any link function. To our knowledge, only [START_REF] Lipovetsky | Analytical closed-form solution for binary logit regression by categorical predictors[END_REF] provides an explicit solution in the special case of a logit regression with categorical predictors. He assumes the response variable follows a Bernoulli distribution with a canonical link function and a particular set-up of intercept.

Then, the paper will continue by the application of such formulas not on classical distributions, but on distributions such as the log-transformed variable has a distribution in the exponential family. Therefore, the choice of probability distributions of this paper is led by two aspects: distributions with positive values and distributions as the log-transformed variable belongs to the exponential family. The considered distributions have heavier tails than the exponential distribution. We choose to study two distributions: the Pareto 1 distribution and the shifted lognormal distribution with fixed threshold parameters. We could have considered log-logistic and GEV distributions being also appropriate in many situations but these distributions do not belong to the exponential family.

Applications of this distribution can be found in various disciplines such as finance, insurance, reliability theory, etc. Here, we are interested with an application to insurance loss modeling. Making insurance tariffs consists in appropriately selecting and transforming explanatory variables so that the prediction fairly estimates the mean of the response variable. When the relation between the transformed response variable and an explanatory variable is not affine, non-linearities is generally accommodated in three ways: binning the variable, adding polynomial terms or using piecewise linear functions in the predictor, see e.g. [START_REF] Goldburd | Generalized linear models for insurance rating[END_REF]. In this paper, we focus mainly on the first approach where continuous variables (typically the age of the policyholder) have been discretized so that explanatory variables are categorical.

More precisely, pricing non-life insurance relies on estimating the claim frequency and the claim severity distributions. The former is generally estimated by a regression model such as Poisson or zero-inflated models. However for modeling claim severity, we commonly split the claim dataset between attritional and atypical claims. A threshold µ is chosen either from the extreme value theory or by expert judgments. A classical GLM such as gamma or inverse-Gaussian is fitted on attritional claim amounts below µ, see e.g. [START_REF] Ohlsson | Non-Life Insurance Pricing with Generalized Linear Models[END_REF]. Atypical claim amounts above µ are not necessarily modeled at all. An empirical rule of the insurance pricing is used to mutualize atypical claims over the portfolio, i.e. the aggregate sum of atypical claims is shared equally among all policies. We aim at providing a regression model for those claims above µ in order to refine this empirical rule.

The threshold µ can also be interpreted in another insurance context. Generally in nonlife insurance, contracts are underwritten with a deductible. This has two consequences: the policyholder will retain the risk of claims below the deductible; and the insurer will only know and be interested in claims above the deductible. In the numerical section, we consider only the example of large claim modeling.

The paper is organized as follows. In Section 2, we present the GLMs. Section 3 provide exact formulas for MLE in the case of categorical explanatory variables. Section 4 is dedicated to the Pareto 1 GLM, while Section 5 is dedicated to the shifted lognormal GLM. Finally, a simulation analysis is provided in Section 6 and an application to an actuarial dataset is carried out in Section 7, before Section 8 conludes.

Preliminaries on Generalized Linear Models

In this section, we consider the estimation problem in GLMs. We consider deterministic exogenous variables x 1 , . . . , x n , with

x i = (x (1) i , . . . , x (p) i ) ∈ R p for i = 1, . . . , n.
In the following, for the sake of clarity, bold notations are reserved for vector of R p and bold notations with an underline are reserved for vector of R n . The index i ∈ I = {1, . . . , n} is reserved for the observations, while the indexes j, k, l are used for the explanatory variables.

In this setting, the sample Y = (Y 1 , . . . , Y n ) is composed of real-valued independent random variables; each one belongs to a family of probability measures of one-parameter exponential type with respective parameters λ 1 , . . . , λ n valued in Λ ⊂ R.

Precisely, the likelihood L associated to the statistical experiment generated by

Y i , i ∈ I verifies log L(ϑ | y i ) = λ i (ϑ)y i -b (λ i (ϑ)) a(φ) + c(y i , φ), y i ∈ Y ⊂ R, (1)
and -∞ if y i / ∈ Y, where a : R → R, b : Λ → R and c : X × R → R are known real-valued measurable functions and φ is the dispersion parameter, e.g. McCullagh and Nelder (1989, Section 2.2).

In Equation ( 1), the parameters λ 1 , . . . , λ n depend on a finite-dimensional parameter ϑ ∈

Θ ⊂ R p . Direct computations lead to b (λ i (ϑ)) = E ϑ Y i and b (λ i (ϑ))a(φ) = Var ϑ Y i .
(2)

Using a twice continuously differentiable and bijective function g from b (Λ) to R, the GLM are defined by assuming the following relation between the expectation E ϑ Y i and the predictor

g(b (λ i (ϑ))) = x i , ϑ = η i , for all ϑ ∈ Θ,
where η i are the linear predictors and ., . denotes the scalar product. In other words, the bijective function = (b ) -1 • g -1 is setted; then we have

λ i (ϑ) = (η i ). ( 3
)
We summarize with the following relations

X × Θ .,. -→ D -1 Λ,
where D is the space of linear predictor and X the possible set of value of x i for i ∈ I. Here is chosen and, consecutively Θ, Λ and D must be set. The parameter ϑ ∈ Θ ⊂ R p is to be estimated and g is called the link function in the regression framework. We talk of canonical link function, when is the identity function.

Let us compute the log-likelihood of y = (y 1 , . . . , y n ):

log L(ϑ | y) = n i=1 y i (η i ) -b ( (η i )) a(φ) + n i=1 c(y i , φ), (4)
with b, h and being respectively defined in (1) and (3). Here, the vector of the parameters ϑ is unknown. If the model is identifiable, it can be shown that the sequence of MLE ( ϑ n ) n≥1 defined by ϑ n = arg max ϑ∈Θ L(ϑ | y) asymptotically exists and is consistent (for example Fahrmeir and Kaufmann, 1985, Theorem 2, 4).

The MLE ϑ n , if it exists, is the solution of the non linear system

S j (ϑ) = 0, j = 1, . . . , p, (5)
with S j (ϑ) are the component of the Score vector defined by

S j (ϑ) = 1 a(φ) n i=1 x (j) i (η i ) (y i -b ( (η i ))) .
It is worth mentioning that for a small data set (small n) or large number of explanatory variables, the existence of the MLE is not guaranteed. Note that the MLE ϑ n does not depend on the value of the dispersion parameter φ. Indeed, the dispersion parameter is estimated in a second step using the sum of square residuals or the log-likelihood, see e.g. (McCullagh and Nelder, 1989, Chap. 9).

In a general setting, the system (5) does not have a closed-form solution and GLMs are generally fitted using a Newton-type method, such that an iteratively re-weighted least square (IWLS) algorithm also refereed to Fisher Scoring algorithm, see e.g. [START_REF] Mccullagh | Generalized linear models[END_REF].

In the case of categorical explanatory variables described later on, the non-asymptotic existence of the MLE depends on the conditional distribution and the chosen link function (see Examples 1, 2 and 3 on Section 4.2).

A closed-form MLE for categorical explanatory variables

In any regression model, categorical or nominal explanatory variables have to be coded since their value is a name or a category. When the possible values are unordered, it is common to use a binary incidence matrix or dummy variables where each row has a single unity in the column of the class to which it belongs. In the case of ordered values, a contrast matrix has to be used, see e.g. [START_REF] Venables | Modern Applied Statistics with S[END_REF].

A single explanatory variable

Let us first consider the case of a single categorical explanatory variable. That is p = 2 and for all i ∈ I, x

i = 1 is the intercept and x

i takes values in a set of d modalities {v 1 , . . . , v d } with d > 2. We define the incidence matrix (x

(2),j i ) i,j where x (2),j i = 1 x (2)
i =v j is the binary dummy of the jth category for i ∈ I and j ∈ J = {1, . . . , d}. From this incidence matrix, we compute the number of appearance m j > 0 of the jth category and y (j) n the mean value of y taking over the jth category by

m j = n i=1 x (2),j i , j ∈ J and y (j) n = 1 m j n i=1 y i x (2),j i , j ∈ J.
By construction, this incidence matrix has rows that sum to 1. Therefore if we use the combination of the incidence matrix with a 1-column for the intercept (x

i , x

(2),j i

) i,j : a redundancy appears. We must choose either to use no intercept, to drop one column for a particular modality of x

(2) i , or to use a zero-sum condition on the parameters. We investigate below these three options in a single framework.

Consider the following GLM for the explanatory variables x

(1)

i , x (2),1 i , . . . , x (2),d i assuming that g(EY i ) = ϑ (1) + d j=1 x (2),j i ϑ (2),j , i ∈ I, (6)
where ϑ = (ϑ (1) , ϑ (2),1 , . . . , ϑ (2),d ) is the unknown vector parameters. The model being not identifiable, we impose exactly one linear equation on ϑ R, ϑ = 0,

with R = (r 0 , r 1 , . . . , r d ) any real vector of size d + 1. A theorem and two corollaries are given below and corresponding proofs are postponed to Appendix A.1.

Theorem 3.1. Suppose that for all i ∈ I, Y i takes values in b (Λ). If the vector R is such that d j=1 r j -r 0 = 0, then there exists a unique, consistent and explicit MLE ϑ n = ( ϑ n,(1) , ϑ n,(2),1 , . . . , ϑ n,(2),d ) of ϑ given by

ϑ n,(1) = d k=1 r k g(Y (k) n ) d k=1 r k -r 0 , ϑ n,(2),j = g(Y (j) n ) - d k=1 r k g(Y (k) n ) d k=1 r k -r 0 , j ∈ J. (8) Note that if Y (j) n does not belong to b (Λ), g(Y (j)
n ) and hence ϑ n,(l),j are not defined.

We give below the three most common examples of linear constraint, some details of these calculus are given in Appendix A.1.

Example 3.1. No-intercept model

The no-intercept model is obtained with R = (1, 0, . . . , 0) leading to ϑ (1) = 0. Therefore the unique, consistent and explicit MLE ϑ n of ϑ is

ϑ n,(1) = 0, ϑ n,(2),j = g Y (j) n , j ∈ J. ( 9
)
Example 3.2. Model without first modality The model without first modality is obtained with R = (0, 1, . . . , 0) leading to ϑ (2),1 = 0. Therefore, the unique, consistent and explicit MLE ϑ n of ϑ is

ϑ n,(1) = g Y (1) n , ϑ n,(2),1 = 0, ϑ n,(2),j = g Y (j) n -ϑ n,1 , j ∈ J \ {1}. Example 3.3. Zero-sum condition
The zero-sum model is obtained with R = (0, 1, . . . , 1) leading to d j=1 ϑ (2),j = 0. Therefore, the unique, consistent and explicit MLE ϑ n of ϑ is

ϑ n,(1) = 1 d d k=1 g Y (k) n , ϑ n,(2),j = g Y (j) n -ϑ n,1 , j ∈ J.
Remark 3.1. In Theorem 3.1, it is worth noting that the value of ϑ n does not depend on the distribution of the Y i . This fact was stated in [START_REF] Goldburd | Generalized linear models for insurance rating[END_REF] but without any proof.

Remark 3.2. The three different parametrizations (Examples 3.1, 3.2 and 3.3) depends on the type of application and on the modeler choice. In statistical software, there is a default choice: for instance in the statistical software R, the model without the first modality is the default parametrization (see functions lm(), glm() by [START_REF] Core | R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing[END_REF]). The first option without intercept may be justified when no group can be chosen as the reference group.

Remark 3.3. When g is the identity function, the first and third options (Examples 3.1 and 3.3) can be interpreted as a generalized analysis of variance (ANOVA) for Y i with respect to groups defined by the explanatory variable x (2) . Even for non-Gaussian random variables, some applications may justify these options.

Remark 3.4. The case when there is no explanatory variable, i.e. Y i are identically distributed, cannot be obtained with Equation (8). But in that case, we get with similar arguments ϑ n,(1) = g(Y n ), which is in line with Example 6.3.12 of [START_REF] Lehmann | Theory of Point Estimation[END_REF] when g is the canonical link, see Appendix A.1.

Remark 3.5. Despite the distribution of Y (j)

n still belongs to the exponential family, the bias of ϑ n,(1) and ϑ n,(2),j cannot be determined for a general link function g. However, we can show the consistency of the estimator and an asymptotic confidence interval, see Appendix A.1. In the following, we will investigate the bias and the distribution of the MLE for some special cases of distributions and link functions.

Theorem 3.1 has two interesting corollaries which give some clues on the choice of the link function g. This corollary tempers the importance of the link function since it will not affect the predicted moments in the case of a single explanatory variable.

Corollary 3.1. The value of the log-likelihood defined in (4) taken on the exact MLE ϑ n (if it exists) given by (8), under constraint (7), does not depend on the link function g. More precisely, we have ∀i ∈ I,

(η i ) = (b ) -1 (y (j) n ) for j ∈ J such that x (2),j i = 1 and log L( ϑ n | y) = 1 a(φ) d j=1 i,x (2),j i =1 y i b y (j) n -b b y (j) n + n i=1 c(y i , φ),
with b = (b ) -1 . Therefore, the criteria AIC and BIC are also independent of the link function g. The estimator of φ is obtained by maximizing log L( ϑ n | y) with respect to φ given a, b, c functions.

Corollary 3.2. The predicted mean and predicted variance for the ith individual is estimated by EY i = b ( ( η i )) and VarY i = a( φ)b ( ( η i )) respectively using (2). Both estimates do not depend on the link function g and the predicted mean does not depend on the function b. More precisely, when v j is the modality of the ith individual (i.e. x

(2),j i = 1), the predicted mean and predicated variance are given by

EY i = y (j) n , VarY i = a( φ)b • (b ) -1 y (j) n .
Corollary 3.2 may be surprising because the predicted mean does not depend on the conditional distribution of Y 1 , . . . , Y n . The predicted mean is just the mean of the response variable taken over the class j i.e. observations y i such that the covariate x

(2) i takes the modality v j . The formula (8) of the MLE ϑ n can be reformulated as

ϑ n = Q R -1 g( Ȳ ) 0 with Q is the d × (1 + d) matrix defined by Q = (A 0 , A 1 )
, with A 0 is the ones vector of size d, A 1 the identity matrix of size d, and g( Ȳ ) the vector (g(

Ȳ (1) n ), . . . , g( Ȳ (d) n )). We have d j=1 r j -r 0 = 0 ⇔ rank Q R = d + 1.
With this formulation, the estimator of ϑ is

ϑ n = (Q Q + R R) -1 Q g( Ȳ )
. This general form is particularly useful in the case of two categorical variables of the next subsection.

Two explanatory variables

Now, we consider the case of two explanatory categorical variables. That is p = 3 and for all i = 1, . . . , n, x

(1) i

= 1 is the intercept and x

(2) i , x

(3) i take values in {v j1 , . . . , v jd j } with d 2 , d 3 modalities respectively. We define by x

(2),k i and x

(3),l i

, k ∈ K = {1, . . . , d 2 } and l ∈ L = {1, . . . , d 3 } the binary dummies of the kth and lth resp. categories, m (j) k > 0 the number of appearance of the kth modality of the jth variable, j = 1, 2, m k,l the number of appearance of the kth and lth category simultaneously and y (k,l) n the mean value of y taking over the kth and lth categories. That is

Dummy Frequency Mean Index x (2),k i = 1 x (2) i =v 2k m (2) k = n i=1 x (2),k i y (2),k n = 1 m (2) k n i=1 y i x (2),k i k ∈ K x (3),l i = 1 x (3) i =v 3l m (3) l = n i=1 x (3),l i y (3),l n = 1 m (3) l n i=1 y i x (3),l i l ∈ L x (k,l) i = x (2),k i x (3),l i m k,l = n i=1 x (k,l) i y (k,l) n = 1 m k,l n i=1 y i x (k,l) i (k, l) ∈ K × L where y (k,l) n is computed over KL = (K × L) \ {(k, l) ∈ K × L; m k,l = 0}. Set d 2,3 = #KL , for l ∈ L, K l = {k ∈ K; m k,l > 0}, d (3),l = #K l and for k ∈ K, L k = {l ∈ L; m k,l > 0}, d (2),k = #L k .
Note that (m k,l ) kl are absolute frequencies of the contingency table resulting from crossclassifying factors and can be computed very easily. Be careful that KL is not equal to

K × L but l∈L K l = k∈K L k = KL , and d 2,3 = d 2 d 3 -r, where r = #{(k, l) ∈ K × L; m k,l = 0}. Let Q be the d 2,3 × (1 + d 2 + d 3 + d 2,3 ) real matrix defined by Q = (A 0 , A 1 , A 2 , A 12 ) with A 0 = 1 d 2,3 the d 2,3 × 1 ones matrix; A 1 = (diag(1 d (2),k )) k∈K , the d 2,3 × K diagonal block matrix of ones vector of size d (2),k ; A 2 = (I ,k d 3 ) k∈K , the d 2,3 × L matrix where I ,k d 3 is the identity matrix of size d 3 without rows l for which m k,l = 0; A 12 = I d 2,3 the d 2,3 × d 2,3 identity matrix.
Consider the following GLM for explanatory variables x

(1) i , x (2),j i , x (3),j i g(EY i ) = ϑ 1 + d 2 k=1 x (2),k i ϑ (2),k + d 3 l=1 x (3),l i ϑ (3),l + (k,l)∈KL x (k,l) i ϑ kl , (10)
where

ϑ (1) , (ϑ (2),k ) k∈K , (ϑ (3),l ) l∈L , (ϑ kl ) (k,l)∈KL are the d 2 + d 3 + d 2,3 + 1 unknown parameters.
Again at this stage, the model is not identifiable because of the redundancy on the vectors (x

(2),k 1 , . . . , x (2),k n ), k ∈ K, the vectors (x (3),l 1 , . . . , x (3)
,l n ), l ∈ L and the ones vector. As previously, we need to impose q ≥ 1 + d 2 + d 3 linear constraints on the vector parameters ϑ

Rϑ = 0 q , ( 11
)
where

R is a q × (1 + d 2 + d 3 + d 2,3 ) real matrix of linear contrasts, with rank(R) = 1 + d 2 + d 3
and 0 q the zeros vector of size q. Again, the proofs of the following theorem and corollaries are postponed to Appendix A.2.

Theorem 3.2. Suppose that for all i ∈ {1, . . . , n}, Y i takes values in b (Λ). Under constraint (11) and if R such that (Q , R ) is of rank d 2,3
, there exists a unique, consistent and explicit MLE ϑ n of ϑ given by

ϑ n = (Q Q + R R) -1 Q g( Ȳ ), (12)
where the vector g( Ȳ ) is

((g( Ȳ (k,l) n )) l∈L k ) k∈K .
Example 3.4. No intercept and no single-variable dummy The model with no intercept and no single-variable dummy is ϑ 1 = 0 and ϑ (2),k = ϑ (3),l = 0 ∀k ∈ K ∀l ∈ L. Therefore, the unique, consistent and explicit MLE ϑ n of ϑ is

ϑ n,kl = g Y (k,l) n , (k, l) ∈ KL .
k ϑ (2),k = l∈L m (3) l ϑ (3),l = 0, ∀l ∈ L, k∈K l m k,l ϑ kl = 0, ∀k ∈ K, l∈L k m k,l ϑ kl = 0.
Therefore, the unique, consistent and explicit MLE ϑ n of ϑ is

                         ϑ n,(1) = 1 n (k,l)∈KL m k,l g Y (k,l) n ϑ n,(2),k = 1 m (2) k l∈L k m k,l g Y (k,l) n -ϑ n,1 , k ∈ K ϑ n,(3),l = 1 m (3) l k∈K l m k,l g Y (k,l) n -ϑ n,1 , l ∈ L ϑ n,kl = g Y (k,l) n -ϑ n,(2),k -ϑ n,(3),l -ϑ n,1 , (k, l) ∈ KL .
Remark 3.6. The MLE of the model with only main effects for two explanatory variables defined as g

(EY i ) = ϑ 1 + d 2 k=1 x (2),k i ϑ (2),k + d 3 l=1 x (3),l i ϑ (3),l
does not present such explicit formula whatever the matrix R of rank 2. In that case, the MLE does not solves a least square problem under a linear constraint, see Appendix A.2. In the special case of logit-regression, Lipovetsky (2015) also notice that least square estimation does not coincide with the MLE.

For simplicity, we consider only the cases of one and two explanatory categorical variables. With a higher number of explanatory variables, we can perform a similar analysis to obtain an explicit solution of the MLE. As for one explanatory variable, Theorem 3.2 has two interesting corollaries on the value of the log-likelihood and the predicted moments.

Corollary 3.3. The value of log-likelihood defined in (4) taken on the exact MLE ϑ n (if it exists) given by (12), under constraint (11), does not depend on the link function g. More precisely, we have ∀i ∈ I,

(η i ) = (b ) -1 (y (k,l) n ) for l ∈ L and k ∈ K such that x (2),j i = 1 and x (3),k i = 1 and log L( ϑ n | y) = 1 a(φ) (k,l)∈KL i∈ Ĩ y i b y (k,l) n -b b y (k,l) n + n i=1 c(y i , φ), with b = (b ) -1 and Ĩ = {i ∈ I, x (2),k i = x (3),l i = 1}.
EY i = y (k,l) n , VarY i = a( φ)b • (b ) -1 (y (k,l) n ).
In the next two sections, we apply previous theorems and corollaries to two particular distributions: Pareto 1 and lognormal distribution. Our results do not only apply to continuous distributions but also for discrete distributions. But we choose these two distributions in order to model insurance losses.

4 GLM for Pareto I distribution with categorical explanatory variables

Characterization

Consider the sample Y = (Y 1 , . . . , Y n ) composed of independent Pareto Type 1. Precisely, we assume that the independent random variables Y 1 , . . . , Y n are Pareto with known threshold parameter µ and respective shape parameter (depending on the unknown parameter ϑ) λ 1 (ϑ), . . . , λ n (ϑ) ∈ Λ = (0, ∞) . The density f of Pareto distribution with threshold and shape parameter µ and λ i (ϑ), i ∈ I is

f (y) = λ i (ϑ) µ λ i (ϑ) y λ i (ϑ)+1 , y ∈ Y = [µ, ∞), (13)
and 0 if y < µ.

We recall that for the Pareto Type 1 distribution

EY i = λ i (ϑ)µ λ i (ϑ) -1 < +∞, iff λ i (ϑ) > 1 and EY 2 i = λ i (ϑ)µ 2 λ i (ϑ) -2 < +∞, iff λ i (ϑ) > 2.
Unlike the known parameter µ, the parameter ϑ is to be estimated. These closed-form formulas are particularly useful in an insurance context since the expectation and the variance are used in most premium computation. For instance, EY is the pure premium and for γ > 0, EY +γVarY is the variance principle (see [START_REF] Bühlmann | A course in credibility theory and its applications[END_REF], Section 1.2.2).

In the following, instead of Y we consider the sample Z = (T (Y 1 ), . . . , T (Y n )). With the re-parametrization z i = T (y i ) = -log (y i /µ), i ∈ I, this distribution belongs to the exponential family as defined in (1), with

a(φ) = 1, b(λ) = -log(λ), and c(z, φ) = 0, z ∈ T (Y) = R -, λ ∈ Λ. ( 14
)
In particular, for the Pareto I distribution, there is no dispersion parameter. It is also worth mentioning that -Z i is exponential with parameter λ i (ϑ). Consecutively, all moments of Z i exist and are given by

E(Z i ) m = (-1) m m!/λ i (ϑ) m , m ∈ N * .
Consider the regression model with a link function g, a response variable Y i Pareto I distributed where

Z i = -log (Y i /µ) and g (EZ i ) = ϑ (1) + x (2),1 i ϑ (2),1 + . . . + x (2),d i ϑ (2),d = x i , ϑ , i ∈ I (15)
with for i ∈ I,

x i = (1, x (2),1 i , . . . , x (2)
,d i

) T are the covariate vectors and ϑ = (ϑ (1) , ϑ (2),1 . . . , ϑ (2),d ) T is the unknown parameter vector.

The choice of the link function g appearing in ( 15) is a crucial point. Let us start with the canonical link. That is, the chosen function g so that = (b ) -1 • g -1 is the identity function. For our Pareto model, g(t) = -1 t since b (λ) = -1 λ . From (4), the choice of the canonical link function imposes constraints on the linear predictor space D ⊂ Λ ⊂ (0, +∞) in that case. Since D results from the scalar product of ϑ parameters with explanatory variables y i , some negative values might be produced when the covariables take negative values. This make the choice of the canonical link inappropriate.

In order to remedy this issue, we can choose a link function such that the values of the function falls in Λ ⊂ (0, ∞). A natural choice is (η) = exp(η) which guarantees a positive parameter. The choice (η) = exp(η)+1 guarantee a finite expectation for the random variables Y i , i = 1, . . . , n. We summarize in Table 1 the tested functions in our application in Section 7, and in Table 2, the spaces given a link function.

(η i ) g -1 (t) g(t) canonical η i - 1 t - 1 t log-inv e η i -e -t log(- 1 t ) shifted log-inv e η i + 1 - 1 e t + 1 log(- 1 t -1)
ϑ × X i space D space Λ b (Λ) unspecified ϑ × X i ⊂ R p × R p .,.
-→

D ⊂ R -1 Λ ⊂ (0, +∞) canonical {(ϑ, x i ) ∈ R p ×R p , ϑ, x i >0}
.,.

-→ (0, +∞)

id id (0, +∞) (-∞, 0) log-inv R p × R p .,. -→ R log exp (0, +∞) (-∞, 0) shifted log-inv R p × R p .,.
-→ R log(x-1) exp(x)+1

(1, +∞) (-1, 0)

Estimation for categorical exogenous variables

Consider the case of one categorical exogenous variable. We expose the case of the re-parametrization without intercept, i.e. R, ϑ = 0 with R = (1, 0, . . . , 0). Let ϑ n the MLE defined in (9), if it exists, of ϑ. Using the following equality of Corollary 3.1

d j=1 i,x (2),j i =1 z i b z (j) n = n, with b = (b ) -1 ,
the log likelihood evaluated on ϑ n for both the transformed sample z and the original sample y with one categorical exogenous variable is

log L( ϑ n | z) = n - d j=1 m j log(-z (j) n ), log L( ϑ n | y) = log L( ϑ n | z) - n i=1 log(y i ). ( 16
)
The log-likelihood computation is detailed on Appendix B.

The second remark is that -Z i are exponentially distributed E( (η i )). Hence for j ∈ J, the estimators ϑ n,(2),j of ϑ (2),j are known transforms of a gamma random variable Ga(m j , m j (ϑ 2,j )).

Below we analyze the choice of the link functions considered in Table 1 in Examples 4.1, 4.2, 4.3 and plotted in Figure 4a. In the special case of canonical Pareto model, because z i < 0 for all i ∈ I, we have z (j) n ∈ b (Λ) = (-∞, 0) for all j ∈ J (g and Λ are respectively defined in Tables 1 and2). With no-intercept using Equation (9), the MLE is

ϑ n,(2),j = -m j n i=1 x (2),j i Z i -1 = - 1 Z (j) n , j ∈ J.
Hence, for j ∈ J, ϑ n,(2),j follows an Inverse Gamma distribution with shape parameter m j and rate parameter m j ϑ (2),j , see e.g. (Johnson et al., 2000, Ch. 17). We can compute the moments of the Inverse Gamma distribution, for m j > 2, E ϑ n,(2),j = m j m j -1 ϑ (2),j , and Var ϑ n,(2),j = m 2 j (m j -1) 2 (m j -2) ϑ 2

(2),j , j ∈ J.

An unbiased estimator of ϑ (2),j is then ϑ n,(2),j = m j -1 m j ϑ n,(2),j which has a lower variance Var ϑ n,(2),j = ϑ 2

(2),j m j -2 ≤ Var ϑ n,(2),j , j ∈ J.

A similar bias is also obtained by [START_REF] Bühlmann | A course in credibility theory and its applications[END_REF] in a credibility context. Of course this unbiased estimator is also applicable for two exogenous variables with the first parametrization of the Theorem 3.2. When some modalities (or couple of modalities) aren't much represented, it can be relevant to use this unbiased estimator.

Example 4.2. log-inverse link In the special case of the log-inv Pareto model, we also have z (j) n ∈ b (Λ) = (-∞, 0) for all j ∈ J (g and Λ are respectively defined in Tables 1 and2). With no-intercept using Equation (9), the MLE is

ϑ n,(2),j = -log 1 -m j n i=1 x (2),j i Z i = -log -Z (j) n , j ∈ J.
Here, for j ∈ J, the distribution of -ϑ n,(2),j is the distribution of the log of the gamma distribution with shape m j and rate m j exp(ϑ (2),j ). We can derive moments of this distribution which should not be confused with the log-gamma distribution studied e.g. in [START_REF] Hogg | Loss distributions[END_REF].

Let L = log(G) when G is gamma distributed with shape parameter a > 0 and rate parameter λ > 0. We have by elementary manipulations the moment generating function of L:

M L (t) = Ee tL = Γ(a + t) Γ(a) λ -t , t > -a,
where Γ denotes the usual gamma function. Therefore by differentiating and evaluating at 0, we deduce that the expectation and the variance of L are EL = ψ(a) -log λ and VarL = ψ (a), where the functions ψ and ψ are the digamma and trigamma function, see e.g. Olver et al. (2010). Getting back to our example, we deduce that E ϑ n,(2),j = ϑ (2),j + log m j -ψ(m j ) and Var ϑ n,(2),j = ψ (m j ), j ∈ J.

From Olver et al. ( 2010), we know that log(m j ) -ψ(m j ) tends to 0 as m j tend to infinity. Hence ϑ n,(2),j is asymptotically unbiased, and an unbiased estimator of ϑ (2),j is ϑ n,(2),j = ϑ n,(2),j -(log(m j ) -ψ(m j )), j ∈ J.

Example 4.3. shifted log-inverse link In the special case of the shifted log-inv Pareto model, z (j) n is not necessarily in b (Λ) = (-1, 0) for all j ∈ J (g and Λ are respectively defined in Tables 1 and2). If there is an index j such as z (j) n ≤ -1, the MLE is not defined and we couldn't use the shifted log-inv link with the same incidence matrix.

Nevertheless, for sufficiently large n, for j such that x

(2) i = v j , Z (j)
n → EZ i almost surely, where EZ i = -1/(exp(ϑ 2,j ) + 1) > -1. Hence for sufficiently large n, the conditions of Theorem 3.1 are satisfied. With no-intercept using Equation (9), the MLE (provided it exists) is

ϑ n,(2),j = log m j -n i=1 x (2),j i Z i -1 = log -1/Z (j) n -1 , j ∈ J.
The expectation of ϑ n,(2),j is complex and should be done numerically. However by the strong law of large numbers and the continuity of the link function, ϑ n,(2),j = -log -1/Z (j)

n -1 converge almost surely to log((exp(ϑ (2),j ) + 1) -1) = ϑ (2),j .

Remark 4.1. In Theorem 3.1, the condition Y i takes values in b (Λ) might seem too restrictive. In fact the condition y (j) n ∈ b (Λ) for all j ∈ J is sufficient to define a vector value ϑ n which maximise the likelihood. But ϑ n fails to be a MLE estimator because the random variable g(Y (j) n ) can to be not defined. Nevertheless, when m j tends to infinity for any j ∈ J, the random variables Y i 's defined by (6) such that y (2),j i = 1 are i.i.d. (not only independent) and the strong law of large numbers implies that Y

Model diagnostic

In this paragraph, we propose residuals adapted at the case of Pareto problem. First note that Y i is Pareto I with shape (η i ) and threshold µ and for the parametrization ( 14) -Z i = log(Y i /µ) ∼ E( (η i )). Let define the residuals

R i = -(η i )Z i , i ∈ I.
Hence R 1 , . . . , R n are i.i.d. and have an exponential distribution E(1). The consistency of the MLE makes it possible to say that the estimated residuals R n,i = -( η i )Z i , i ∈ I, with η i = x i , ϑ n are asymptotically i.i.d..

It is also possible to verify the assumption of the Pareto distribution for Y i conditionally to y i with graphical diagnostic as an exponential Quantile-Quantile plot on the residuals R n,i .

In the case of a single explanatory variable, for i ∈ I, the residuals R n,i do not depend on function. Their explicit forms are

R n,i = Z i Z(j) n j such that x (2) i = v j , i ∈ I. (17)
Furthermore, the summation of R n,1 , . . . , R n,n has the surprising property to be deterministic and is exactly equal to n.

5 GLM for shifted log-normal distribution with categorical explanatory variables

Characterization

Secondly, consider the sample Y = (Y 1 , . . . , Y n ) composed of independent shifted log-normal variables respectively with mean λ 1 (ϑ), . . . , λ n (ϑ), dispersion φ = σ 2 and a known threshold µ. The shifted log-normal is also known as the 3-parameter log-normal. Precisely, the density of Y i is

f (y) = 1 (y -µ) √ 2πφ exp - (log(y -µ) -λ i (ϑ)) 2 2φ , y ∈ Y = (µ, ∞), (18)
and 0 for y ≤ µ. It is well known that the lognormal distribution has finite moment, see e.g. [START_REF] Johnson | Continuous Univariate Distributions[END_REF]. In particular, the expectation and the variance are given by

EY i = µ + exp(λ i (ϑ) + φ/2), VarY i = (exp(φ) -1) exp(2λ i (ϑ) + φ).
The transformed sample Z = T (Y ) = (T (Y 1 ), . . . , T (Y n )) with T (y) = log(y -µ) is belongs the exponential family with

a(φ) = φ, b(λ) = λ 2 /2, c(z, φ) = - 1 2 z 2 φ + log(2πφ) , z ∈ R+, λ ∈ R.
It is worth mentioning that Z i are normally distributed with mean λ i (θ) and variance φ. As a consequence, all moments of Z i exists and E(Z i -λ i ) m = (2m)!φ m /(2 m m!) for m even and 0 for m odd. Consider the regression model with a link function g, a response variable Y i lognormally distributed where

Z i = log (Y i -µ) and g (EZ i ) = ϑ (1) + ϑ (2),1 x (2),1 i + . . . + ϑ (2),d x (2),d i = x i , ϑ , i ∈ I (19)
with for i ∈ I,

x i = (1, x (2),1 i , . . . , x (2),d i
) T are the covariate vectors and ϑ = (ϑ (1) , ϑ (2),1 , . . . , ϑ (2),d ) T ∈ R d is the unknown parameter vector.

The choice of the link function g for Equation ( 19) is less restrictive than for the Pareto case. Any differentiable invertible function from R to R will work. Since b (x) = x, the canonical link function is obtained by choosing g such that = id • g -1 = g -1 is the identity function. In other words, the canonical link function is the identity function.

Unlike the previous section, any moment of Y i exist and there is no particular link needed to guarantee their existence. In the numerical application, we will also consider another link function: a real version of the logarithm.

Estimation for categorical exogenous variables

Again we consider the case of categorical variables and without intercept model, that is with a predictor

η i = x (2),1 i ϑ (2),1 + • • • + x (2),d i ϑ (2),d .
In the case of the lognormal dispersion, there is a dispersion to be estimated. The log-likelihood is given by log

L(ϑ | z) = - 1 2φ n i=1 (z i -λ i (ϑ)) 2 - n log(2πφ) 2 .
The components of the MLE of ϑ are given by ϑ n,(2),j = g(z

(j)
n ), j ∈ J, and the estimated log likelihood for the transformed sample z is

log L( ϑ n | z) = - 1 2φ j∈J i,x (2),j i =1 z i -z(j) n 2 - n log(2πφ) 2 .
Maximizing over φ the log-likelihood leads to the empirical variance

φ = 1 n j∈J i,x (2),j i =1 z i -z(j) n 2 . ( 20
)
Hence the estimator of φ is the intra-class variance. This closed-form estimate φ differs from what classical statistical softwares carry out, where the dispersion parameter is estimated by a quasi-likelihood approach.

Using the previous equations, we compute the log likelihood evaluated on φ and on ϑ n for both the transformed sample z and the original sample y with one categorical exogenous variable (Corollary 3.1) is

log L( ϑ n | z) = - n 2 (1 + log(2π φ)), log L( ϑ n | y) = - n 2 (1 + log(2π φ)) - n i=1 z i . ( 21
)
The log-likelihood computation is detailed on Appendix B. Below we analyze the choice of the link functions considered in Table 3 in Examples 5.1, 5.2 and plotted in Figure 4b and with parameter spaces given in Table 4.

Names (η i ) g -1 (t) g(t) canonical η i t t sym. log e η i 1 η i ≥0 + (2 -e -η i )1 η i <0 e t 1 t≥0 + (2 -e -t )1 t<0 log(t)1 t≥1 -log(2 -t)1 t<1
ϑ × X i space D space Λ canonical R p × R p .,. -→ R id id R sym. log R p × R p .,. -→ R l -1 g lg R
Example 5.1. canonical link With the canonical link function, there is no issue between the parameter space and the linear predictor space since D = Λ = R. With no-intercept using Equation (9), the MLE is

ϑ n,(2),j = 1 m j n i=1 x (2,j) i Z i = Z (j) n , j ∈ J.
Hence, the distribution ϑ n,(2),j is simply a normal distribution with mean ϑ (2),j and variance φ/m j . Therefore, the MLE is unbiased and converges in almost surely to ϑ (2),j .

Example 5.2. symmetrical log link We consider a central symmetry of the logarithm function given in Table 3 leading to l g (x) = e x 1 x≥0 + (2 -e -x )1 x<0 . With this symmetrical log link function, there is no issue between the parameter space and the linear predictor space since again D = Λ = R. With no-intercept using Equation (9), the MLE is

ϑ n,(2),j = log Z (j) n 1 Z (j) n ≥1 -log 2 -Z (j) n 1 Z (j) n <1
, j ∈ J. The expectation of ϑ n,(2),j is complex and should be done numerically. However by the strong law of large numbers and the continuity of the link function, ϑ n,(2),j = l g (Z

(j) n ) converge almost surely to l g (EZ (j) n ) = ϑ (2),j .

Model diagnostic

In this paragraph, we give some details about residuals in the lognormal case. As already said, the transformed variables Z i = log(Y i -µ) is normally distributed with mean (η i ) and variance φ. Let define the residuals

R i = Z i -(η i ) √ φ , i ∈ I.
Hence R 1 , . . . , R n are i.i.d. and have a normal distribution N (0, 1). The consistency of the MLE makes it possible to say that the R n,i = Z i -( η i ) √ φ , i ∈ I, with η i = x i , ϑ n are asymptotically i.i.d. Furthermore, the summation of R n,1 , . . . , R n,n is exactly equal to 0.

It is also possible to verify the assumption of the lognormal distribution for Y i conditionally to x i with graphical diagnostic as a normal Quantile-Quantile plot on the residuals R n,i .

Simulations study

This section is devoted to the simulation study: all computations are carried out thanks to the R statistical software [START_REF] Core | R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing[END_REF]. The first part of the simulation analysis consists in assessing the uncertainty of the MLE with the two different approaches: either the explicit formula given in Example 4.1 or the IWLS algorithm described in [START_REF] Mccullagh | Generalized linear models[END_REF].

Therefore, the simulation process has the following steps. Firstly, given the parameter number p, we simulate n random variables which are Pareto 1 distributed in the case of no intercept and a canonical link (see Example 4.1). Secondly, we check that the fitted coefficients by both methods have identical values. Thirdly, we compute the exact confidence interval using the result of Example 4.1 and the asymptotic MLE confidence interval resulted from the IWLS algorithm.

Figure 1 shows the estimated parameters when p = 3 for different values of n. That is we plot ϑ n,(2),1 , ϑ n,(2),2 , ϑ n,(2),3 as solid lines for n = 100, 300, 500, 700, 900, 1000, 3000, 5000, 7000, 9000, 10000 against the true values ϑ = (2, 3, 4) (horizontal dot-dashed lines) and µ = 150. The confidence intervals are plotted in dashed lines (theoretical) and dotted lines (asymptotic). We observe that for any sample size n, explicit and IWLS produce the same value but the explicit confidence interval is much thinner in the explicit case.

Figure 1: Estimated coefficient ϑ n,(2),i (solid lines) with theoretical confidence intervals (dashed lines) and asymptotic confidence intervals (dotted lines)

The second part of the simulation analysis consists in assessing the gain in terms of computations between the two different approaches: either the explicit formula given in Example 4.1 or the IWLS algorithm described in [START_REF] Mccullagh | Generalized linear models[END_REF].

In Table 5, we provide a complexity analysis of the two procedure : the IWLS algorithm and the explicit solution. Again we simulate with Pareto 1 distributed random variables for a specific number of parameters p, and a known threshold µ = 150 given a sample size n. Then we compute the floating point operation numbers given the size of the input dataset. Two different explanatory variables have been tested so that the parameter number is different: p = 5 or 7. We observe that the explicit solution is far less computer intensive (4000 times faster) than the IWLS algorithm which takes 5 or 6 iterations to reach the solution. Thus having an explicit solution can lead to substantial benefits for very large datasets.

IWLS

Application to large claim modeling

This section is devoted to the numerical illustration on a real dataset (again computations are carried out thanks to the R statistical software R Core Team ( 2019)). In our application, we focus on modeling non-life insurance losses (claim amount) of corporate business lines. Our data set comes from an anonymous private insurer: for privacy reason, amounts have been randomly scaled, dates randomly rearranged, variable modalities renamed. The data set consists of 211,739 claims which occurred between 2000 and 2010. In addition to the claim amount level, various explanatory variables are available. We provide in Table 10 in Appendix C a short descriptive analysis of the two most important variables (risk class and guarantee type with respectively 5 and 7 modalities). Due to the very high value of skewness and kurtosis, we observe that claim amount is particularly heavy tailed.

In the sequel, we consider only large claims which are in our context claims above µ = 340, 000 (in euros). The threshold value has been chosen by expert opinion of practitioners. We refer to e.g. [START_REF] Reiss | Statistical Analysis of Extreme Values[END_REF] for advanced selection methods based on extreme value theory.

A single explanatory variable

Firstly, we consider both Pareto 1 GLM and Shifted log-normal GLM with only one explanatory variable: the guarantee type. We choose Guarantee 1 as the reference level implying that ϑ (2),1 = 0. So, ϑ (1) representing the effect of the reference category and (ϑ (2),j ) j representing the differential effect of categories j relative to the reference category will be estimated through ( 15) and ( 19). Observations y 1 , . . . , y n are observed claim amounts either from Pareto 1 (13) or shifted log-normal (18).

For theses two models, we have many possible choices for the link function g. Naturally, we choose link functions appearing in Tables 1 and3 respectively. In accordance to Corollaries 3.1, the choice of g does not impact values respectively given on ( 16) and ( 21) of the loglikelihoods applied on the MLE of ϑ.

Furthermore, for Pareto GLM, the choice of shifted log-inverse link function seems attractive because it guarantees the existence of EY i . Nevertheless, alternative link functions (canonical or log-inv) allow to construct an unbiased estimator (see Section 4). For shifted log-normal model, the choice of canonical link function is more attractive because it leads to an unbiased and simpler MLE estimator (see Section 5).

Coefficients are estimated by explicit formulas given in Sections 4 and 5. In Table 6, the estimated coefficients are given in the five considered situations. Positive values of ϑ (2),j in the Pareto GLM increase the shape parameter of the Pareto 1 distribution leading to a decrease in heavy-tailedness. Regarding the log-normal model, positive values of ϑ (2),j increase the scale parameter of the log-normal distribution leading to a shrink of the distribution.

Irrespectively of the considered link function, the sign of the fitted coefficients are same except for intercept (Table 6) given a distribution. This convinces us that different model assumptions (i.e. link) do not lead to opposite conclusions on the claim severity. Furthermore from Table 10, we retrieve the fact that all guarantees except Guarantee 2 have heavier tails than the reference Guarantee 1. Whatever the considered link function g, the residuals defined in Section 4.3 by R n,i = -( η i )Z i , i ∈ I, do not depend on and are given by Equation ( 17). We show on Figure 3 (left) the quantile/quantile plots of residuals described on Section 4.3 against the standard exponential distribution. On Figure 2, we observe that the assumption of Pareto 1 for Y 1 , . . . , Y n is better than the log-normal distribution at least for small quantiles. Moreover, comparing the value of the log-likelihood in Table 6, Pareto 1 distribution is also the best choice. In the following, we focus only on the Pareto 1 distribution.

For all coefficient, let us compute the p-values of statistical tests with the null hypothesis ϑ (1) = 0 (Interecept null), and for j ∈ J \ {1}, ϑ (2),j = 0 (no differential effect of the jth Guarantee). Table 7 reports the value of the coefficient, its standard error, the student statistics and the associated p-value. We observe that some modalities of the guarantee variable are

Two explanatory variables

Secondly, we consider the Pareto GLM models and Shifted log-normal GLM models with the two explanatory variables (guarantee and risk class) without intercept nor single-variable (c.f. model ( 10) and example 3.4), that are

g(EZ i ) = (k,l)∈KL ϑ kl x (k,l) i , i ∈ I (22)
with Z i = -log(Y i /µ) for the Pareto 1 modeling Z i = log(Y i -µ) for the shifted log normal modeling and where for (k, l) ∈ KL the unknown parameters ϑ kl represent the effect of the couple of the modalities k and l for the first and the second variable. In theses examples, as it describes in Table 8, we have K = {1, . . . , 7}, L = {1, . . . , 5} but KL = {1, . . . , 7} × {1, . . . , 5} \ {(1, 2), (6, 5)}.

Consider the estimation procedures in ( 22). We compute the claim numbers according Guarantee and Risk in Table 8. This claim number per class might be too short to ensure the existence of the MLE with the shifted log-inv link. We arbitrary choose the simple case of the canonical link and an unbiased estimator is relevant in this context.

The coefficients of the model are estimated using the exact method described in Section 3 and then unbiased in the same way of Example 1. The fitted coefficients are not shown but are available upon request to the authors. Furthermore, we compute the p-values of the statistical test ϑ kl = 0 in Table 9. We observe that most computed p-values are small: either less than 10 -6 or less than 1%. Only 5 on the 33 p-values are above the usual 5% level, corresponding to the couples Guarantee/Risk class (1,5), (2,2), (2,3), (2,5) and(7,5) (claim number of 1,2 or 3). In the two variables setting, the Pareto 1 GLM is thus still relevant. < 10 -6 < 10 -6 < 10 -6 < 10 -6 0.00001 Guarantee 5

< 10 -6 0.00002 < 10 -6 < 10 -6 0.01431 Guarantee 6

< 10 -6 0.00815 0.04550 0.00091 -Guarantee 7

< 10 -6 0.00270 < 10 -6 < 10 -6 0.08326

Conclusion

In this paper, we deal with regression models where the response variable belongs to the general formulation of the exponential family, the so-called GLM. We focus on the estimation of parameters of GLMs and derive explicit formulas for MLE in the case of categorical explanatory variables. In this case, the closed-form estimators do not require any use of numerical algorithms, in particular the well-known IWLS algorithm. This is logical, because in the special setting of categorical variables, a regression model is equivalent to fitting the same distribution on subgroups defined with respect to explanatory variables. Hence, we get back to usual explicit solutions for the exponential family in the i.i.d. case.

Yet we work with one or two explanatory variables for the two derived theorems, the approach can be extended to d categorical variables as long as we consider interactions terms and a zero-sum condition. If we consider main effects only for d categorical variables, the MLE cannot be reformulated as a least-square problem. Nevertheless, having an explicit formula make a clear advantage compared to the IWLS algorithm, particularly for large scale datasets.

The explicit formulas are examplified on two particular positive distributions particularly useful in an insurance context: the Pareto 1 distribution and the shifted log-normal distribution. In both cases, we present typical link functions and derive in most cases the distribution of the MLE. In relevant cases, we also give an unbiased estimator. In the general setting, the exact standard error computation is not available, yet the Delta Method can be used to obtain an asymptotic standard error. Finally, we illustrate the estimation process for both distributions on simulated datasets and an actuarial data set.

For future research, a natural extension is to propose regression models for distribution outside the exponential family. Typically, we could consider the Pareto 1 distribution with unknown threshold and shape parameters. We could also consider generalized Pareto distribution based on the peak over thresholds approach. A natural extension could also be to jointly estimate the shape and the dispersion parameters of the distribution.

The system S(ϑ) = 0 is

         n i=1 (η i ) (y i -b • (η i )) = 0 n i=1 x (2),j i (η i ) (y i -b • (η i )) = 0, ∀j ∈ J. that is            j∈J (ϑ (1) + ϑ (2),j ) n i=1 x (2),j i y i -m j b • (ϑ (1) + ϑ (2),j ) = 0 (ϑ (1) + ϑ (2),j ) n i=1 x (2),j i y i -m j b • (ϑ (1) + ϑ (2),j ) = 0, ∀j ∈ J.
The first equation in the previous system is redundancy, and

S(ϑ) = 0 ⇔ (ϑ (1) + ϑ (2),j ) n i=1 x (2),j i y i -m j b • (ϑ (1) + ϑ (2),j ) = 0, ∀j ∈ J.
Hence if Y i takes values in Y ⊂ b (Λ), and injective, we have

ϑ (1) + ϑ (j) = g(Y (j) n ) ∀j ∈ J. The system (23) is Qϑ = g( Ȳ ) Rϑ = 0. ⇔ Q R ϑ = g( Ȳ ) 0 . (24)
Let us compute the determinant of the matrix M d = Q R . Consider R = (r 0 , r 1 , . . . , r d ). We have

M d = 1 d I d r 0 r =        1 1 0 . . . 1 0 1 0 . . . . . . . . . . . . . . . 1 0 . . . 0 1 r 0 r 1 . . . r d        , with r = r 1 . . . r d , 1 d =    1 . . . 1    .
The determinant can be computed recursively

det(M d ) = r d 1 1 0 . . . 1 0 . . . 0 . . . . . . . . . 1 1 0 . . . 0 - 1 1 0 . . . 1 0 . . . 0 . . . . . . . . . 1 r 0 r 1 . . . r d-1 = (-1) d+1 r d -det(M d-1 ). Since det(M 1 ) = -r 0 + r 1 and det(M 2 ) = -r 2 -(-r 0 + r 1 ) = r 0 -r 1 -r 2 , we get det(M d ) = (-1) d r 0 + (-1) d+1 (r 1 + • • • + r d ) = (-1) d (r 0 -r 1 -• • • -r d ).
This determinant is non zero as long as r 0 = d j=1 r j . Now we compute the inverse of matrix M d by a direct inversion.

1 d I d r 0 r a b C d = I d 0 0 1 ⇔        1 d a + I d C = I d b1 d + I d d = 0 r 0 a + rC = 0 br 0 + rd = 1 ⇔          C = I d - 1 -r 0 +r1 d 1 d r d = 1 -r 0 +r1 d 1 d a = r -r 0 +r1 d b = -1 -r 0 +r1 d Let us check the inverse of M d 1 d I d r 0 r r -r 0 +r1 d -1 -r 0 +r1 d I d -1 d r -r 0 +r1 d 1 d -r 0 +r1 d = 1 d r -r 0 +r1 d + I d -1 d r -r 0 +r1 d -1 d -r 0 +r1 d + 1 d -r 0 +r1 d r 0 r -r 0 +r1 d + r -r1 d r -r 0 +r1 d -r 0 -r 0 +r1 d + r1 d -r 0 +r1 d = I d 0 0 1 . So as long as r 0 = d j=1 r j ϑ n = r -r 0 +r1 d -1 -r 0 +r1 d I d -1 d r -r 0 +r1 d 1 d -r 0 +r1 d g( Ȳ ) 0 = rg( Ȳ ) -r 0 +r1 d g( Ȳ ) -1 d rg( Ȳ ) -r 0 +r1 d .
In an other way, the system ( 24) is equivalent to

(Q , R ) Q R ϑ = Q g( Ȳ ),
and for (Q R) of full rank, the matrix (Q Q+R R) is invertible and ϑ = (Q Q+R R) -1 Q g( Ȳ ).

Examples -Choice of the contrast vector R 1. Taking r 0 = 1, r = 0 leads to -r 0 + r1 d = -1 ⇒ ϑ n = 0 g( Ȳ )

.

2. Taking r 0 = 0, r = (1, 0) leads to

-r 0 + r1 d = 1 ⇒ ϑ n =        g( Ȳ (1) n ) 0 g( Ȳ (2) n ) -g( Ȳ (1) n ) . . . g( Ȳ (d) n ) -g( Ȳ (1) n ))        .
3. Taking r 0 = 0, r = 1 leads to

-r 0 + r1 d = d ⇒ ϑ n =      g( Ȳ ) g( Ȳ (1) n ) -g( Ȳ ) . . . g( Ȳ (d) n ) -g( Ȳ )      , with g( Ȳ ) = 1 d d j=1 g(Y (j) n ).
Proof of Remark 3.4. We have to solve the system

S(ϑ) = 0 ⇔ n i=1 (η) (y i -b • (η)) = 0.
If is injective, the system simplifies to -→ n→+∞ g(g -1 (ϑ (1) + ϑ (j) )) = ϑ (1) + ϑ (j) .

By the Central Limit Theorem (i.e. Y In fact, we must be verified than ( η i ) does not depend on g function. If we consider ϑ n defined by (8), we have Q ϑ n = g(ȳ) , since ϑ n is solution of the system (23), i.e. Q(Q Q+R R) -1 Q = I Using η i = (Q ϑ n ) j for i such that x (2),j i = 1 we obtain ). In the same way, taking summation over K l for l ∈ L and over L k for k ∈ K, we found ϑ (2),k and ϑ (3),l , and then ϑ kl .

( η i ) = d j=1 • g(ȳ (j) n )x (2),j i = d j=1 • -1 • (b ) -1 (ȳ (j) n )x (2),j i = d j=1 (b ) -1 (ȳ (j) n )x (
With main effect only, the system S(ϑ) = 0 is Even with a canonical link (x) = x so that (x) = 1, this system is not a least-square problem for a nonlinear g function.

                   n i=1 (η i )y i = n i=1 g -1 (η i ) (η i ) n i=1 x (3),l i (η i )y i = n i=1 x (3),l i g -1 (η i ) (η i ) ∀l ∈ L n i=1 x (2),k i (η i )y i = n i=1 x (2),k i g -1 (η i ) (η i ),
B Calculus of the Log-likelihoods appearing in Sections 4 and 5

Consider the Pareto GLM described on ( 13) and ( 15). The b function is b(λ) = -log(λ), using corollary 3.1 we have (η i ) = (b ) -1 (z (j) n ) = -(z (j) n ) -1 for j such that x

  Example 3.5. Zero-sum conditions The model with zero-sum conditions assumes k∈K m

  The estimator of φ is obtained by maximizing log L( ϑ n | y) with respect to φ given a, b, c functions.Corollary 3.4. The predicted mean and predicted variance for the ith individual is estimated by EY i = b ( ( η i )) and VarY i = a( φ)b ( ( η i )) respectively using (2). Both estimates do not depend on the link function g and the predicted mean does not depend on the function b. Let v 2k and v 3l be the modalities of the ith individual of the two explanatory variables, i.e. x The predicted mean and variance are given by

i

  converges almost surely to EY i = b ( (η i )) ∈ b (Λ). Hence, the probability P (Y (j) n / ∈ b (Λ)) tends to zero which guarantees the asymptotically existence of the MLE estimator.

Figure 2 :

 2 Figure 2: Quantile-quantile plots of the residuals defined in Section 4.3: (left) for Pareto 1, (right) for shifted log-normal distributions

Figure 3 :

 3 Figure 3: Quantile-quantile plots of the residuals defined in Section 4.3: (left) for one explanatory variable, (right) for two explanatory variables

y

  i -nb • (b ) -1 • g -1 (η) = 0 ⇔ η = g y n ⇔ θ = g y n .Proof of Remark 3.5. Let Y i from the exponential family F exp(a, b, c, λ, φ). It is well known, that the moment generating function of Y i isEe tY i = exp b(λ + ta(φ)) -b(λ) a(φ) .Hence, the moment generating function of the average Y m isM Y m (t) = exp b(λ + t m a(φ)) -b(λ) a(φ) m = exp b(λ + ta(φ)/m) -b(λ) a(φ)/m .So we get back to a known result that Y m belongs to the exponential family F exp (x → a(x)/m, b, c, λ, φ), e.g.[START_REF] Mccullagh | Generalized linear models[END_REF].In our setting, random variables in the average Y (j) n are i.i.d. with functions a, b, c and parameters λ = (ϑ (1) + ϑ (j) ) and φ. And Y (j)n also belongs to the exponential family with the same parameter but with the function ā : x → a(x)/m j . In particular,EY (j) n = b ( (ϑ (1) + ϑ (j) )) = g -1 (ϑ (1) + ϑ (j) ), VarY (j) n = a(φ) m j b ( (ϑ (1) + ϑ (j) )).But the computation of Eg(Y (j) n ) remains difficult unless g is a linear function. By the strong law of large numbers, as m j → +∞, the estimator is consistent since Y

n

  converges in distribution to a normal distribution) and using the Delta Method, we obtain that the following √ m a(φ)b ( (ϑ (1) + ϑ (j) ))g (g -1 (ϑ (1) + ϑ (j) )) 2 .Proof of Corollaries 3.1. The log likelihood of ϑ n is defined bylog L( ϑ n | y) = 1 a(φ) n i=1 (y i ( η i ) -b( ( η i ))) + n i=1 c(y i , φ).

,

  In the same way,EY i = b ( ( η i )) = VarY i = a(φ)b ( ( η i )) = a(φ) d j=1 b • (b ) -1 (ȳ (j) n )x(2),j i .2. Under linear contrasts ( CΣ ), the systemϑ (1) + ϑ (2),k + ϑ (3),l + ϑ kl = g( Ȳ (k,l) n ) ∀(k, l) ∈ KL implies that (k,l)∈KL m k,l (ϑ (1) + ϑ (2),k + ϑ (3),l + ϑ kl ) = (k,l)∈KLm k,l g( Ȳ (k,l)

  ∀k ∈ KThere are 1 + d 2 + d 3 equations for 1 + d 2 + d 3 parameters, but each explanatory variable are colinear. So, the two additional constraints Rϑ = 0 ensures that a solution exist for the remainingd 2 + d 3 -1 parameters. Using k x (2),k i = 1, the second set of equations becomes ∀l ∈ L k∈K (ϑ (1) + ϑ (2),k + ϑ (3),l )ȳ (k,l) n m k,l = k∈K g -1 (ϑ (1) + ϑ (2),k + ϑ (3),l ) (ϑ (1) + ϑ (2),k + ϑ (3),l )m k,lSimilarly, the third set of equations becomes ∀k ∈ Kl∈L (ϑ (1) + ϑ (2),k + ϑ (3),l )ȳ (k,l) n m k,l = l∈L g -1 (ϑ (1) + ϑ (2),k + ϑ (3),l ) (ϑ (1) + ϑ (2),k + ϑ (3),l )m k,l

z

  i /z (j) n -log -z (j) n = n -d j=1 m j log -z (j) n .

Table 1 :

 1 Table of typical link functions for Pareto I

	Names

Table 2 :

 2 Summary of spaces for Pareto I

		Parameter-covariable	Linear predictor	Parameter
	Link name	space

Table 3 :

 3 Table of typical link functions for lognormal

Table 4 :

 4 Summary of spaces for lognormal

		Parameter-covariable	Linear predictor	Parameter
	Link name	space

Table 5 :

 5 Floating point operation number given the size of the dataset

		algorithm		exact method
	n p iter. nb. flop x1000	n p iter. nb. flop x1000 rel. gain
	200 5	6	1206.2	200 5	1	1.6	753.9
	400 5	5	4010.4	400 5	1	3.2	1253.2
	600 5	6	10818.6	600 5	1	4.8	2253.9
	800 5	5	16020.8	800 5	1	6.4	2503.2
	1000 5	5	25026	1000 5	1	8	3128.2
	1200 5	5	36031.2	1200 5	1	9.6	3753.2
	1400 5	5	49036.4	1400 5	1	11.2	4378.2
	1600 5	5	64041.6	1600 5	1	12.8	5003.2
	1800 5	5	81046.8	1800 5	1	14.4	5628.2
	200 7	6	1686.2	200 7	1	2	843.1
	400 7	6	6732.4	400 7	1	4	1683.1
	600 7	6	15138.6	600 7	1	6	2523.1
	800 7	6	26904.8	800 7	1	8	3363.1
	1000 7	6	42031	1000 7	1	10	4203.1
	1200 7	6	60517.2	1200 7	1	12	5043.1
	1400 7	6	82363.4	1400 7	1	14	5883.1
	1600 7	6	107569.6 1600 7	1	16	6723.1
	1800 7	6	136135.8 1800 7	1	18	7563.1

Table 6 :

 6 Coefficients for the guarantee variable

	Model		Pareto 1		Shifted log normal
	Variable	canonical	loginv	shifted.loginv	canonical	symlog
	Intercept	1.89	0.64	-0.11	11.75	2.46
	Guarantee 2	0.04	0.02	0.04	0.10	0.01
	Guarantee 3	-0.67	-0.43	-1.36	0.75	0.06
	Guarantee 4	-0.86	-0.60	-3.13	1.04	0.08
	Guarantee 5	-0.71	-0.47	-1.55	0.72	0.06
	Guarantee 6	-0.42	-0.25	-0.63	0.42	0.04
	Guarantee 7	-0.48	-0.29	-0.78	0.59	0.05
	log likelihood -14507.53 -14507.53	-14507.53	-14517.37	-14517.37

Table 7 :

 7 Statistics and p.values for the tests ϑ (1) = 0 and ϑ j = 0, j ∈ J in the Pareto GLM model for the log-inverse link.

		Estimate Std. Error z value Pr(>|z|)
	Intercept	0.6391	0.1644 3.8877	0.0001
	Guarantee 2	0.0214	0.2219 0.0966	0.9230
	Guarantee 3	-0.4332	0.1950 -2.2217	0.0263
	Guarantee 4	-0.6009	0.1708 -3.5180	0.0004
	Guarantee 5	-0.4660	0.1805 -2.5817	0.0098
	Guarantee 6	-0.2485	0.2295 -1.0827	0.2789
	Guarantee 7	-0.2949	0.1834 -1.6075	0.1080

Table 8 :

 8 Number of claim per Guarantee and per Risk class.

	Claim number Risk class 1 Risk class 2 Risk class 3 Risk class 4 Risk class 5
	Guarantee 1	39	0	4	6	1
	Guarantee 2	26	2	3	16	3
	Guarantee 3	48	7	11	29	4
	Guarantee 4	232	40	75	147	20
	Guarantee 5	68	18	36	72	6
	Guarantee 6	24	7	4	11	0
	Guarantee 7	94	9	22	57	3

Table 9 :

 9

	Guarantee 1	< 10 -6	-	0.04550	0.01431	0.31731
	Guarantee 2	< 10 -6	0.15730	0.08301	0.00006	0.08326
	Guarantee 3	< 10 -6	0.00815	0.00091	< 10 -6	0.04550
	Guarantee 4

p-values for the tests ϑ kl = 0, (k, l) ∈ KL in (22) for the canonical link.

p-values Risk class 1 Risk class 2 Risk class 3 Risk class 4 Risk class 5

Acknowledgments

The authors thank Vanessa Desert for her active support during the writing of this paper. The authors are also very grateful for the useful suggestions of the two referees. This work is supported by the research project "PANORisk" and Région Pays de la Loire (France). This is a post-peer-review, pre-copyedit version of an article published in Computational Statistics. The final authenticated version is available online at: http://dx.doi.org/10. 1007/s00180-019-00918-7.

A Proofs of Section 3

A.1 Proof for the one-variable case Proof of Theorem 3.1. We have to solve the system S(ϑ) = 0 Rϑ = 0.

(23)

A.2 Proof for the two-variable case

Proof of Theorem 3.2. The system S(ϑ) = 0 is

that is

The system have exactly 1 + d 2 + d 3 redundancies, and S(ϑ) = 0 reduces to

Hence the system has rank KL and if Y i takes values in Y ⊂ b (Λ), and injective, we have

In the same way of proof of Theorem 3.1, we have to solve

that is, because QQ + RR is full rank, in the same way of proof of Theorem 3.1

In that case, the MLE solves a least square problem with response variable g( Ȳ ), explanatory variable Q under a linear constraint R.

1. Under linear contrasts ( C0 ), the model ( 10) is equivalent to model ( 6) with J = KL modalities. Hence the solution is evident.

Compute the original log likelihood of Pareto 1 distribution:

Hence with

Now consider the shifted log-normal GLM described on ( 18) and ( 19). Here, the b function is b(λ) = λ 2 /2, hence using Corollary 3.1, we have

n for j such that x

(2),j i = 1 and equation ( 21) holds. Let us compute the original log likelihood of the shifted log normal distribution: