# Talk at the Economics Seminar at UQAM

Wednesday, at lunch time, I will give a talk on insurance pricing at the (internal) economics seminar, at UQAM. Slides are now online.

In May, I will also give a seminar on the same topic at the Computational Science seminar…

# Annual UCSB InsurTech Summit

Just a quick post to mention that an Insurtech Summit will be organized in May 2019, on Friday 3rd, by Mike Ludkovski, and I will be there, with Francois Millard (Vitality Group), Adam Tashman (Carpe Data, Santa Barbara), Emil Valdez (University of Connecticut), and Howard Zail (Elucidor, LLC, New York City). That will be nice… I will actually also give a talk on the Monday before at the actuarial seminar !

# Foundations of Machine Learning, part 5

This post is the nineth (and probably last) one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 8 is online here.

## Optimization and algorithmic aspects

In econometrics, (numerical) optimization became omnipresent as soon as we left the Gaussian model. We briefly mentioned it in the section on the exponential family, and the use of the Fisher score (gradient descent) to solve the first order condition $\mathbf{X}^T W(\beta)^{-1})[y-\widehat{y}]=\mathbf{0}$. In learning, optimization is the central tool. And it is necessary to have effective optimization algorithms, to solve problems (described previously) of the form: $$\widehat{\beta}\in\underset{\beta\in\mathbb{R}^p}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,\beta_0+\mathbf{x}^T\beta)+\lambda\Vert\boldsymbol{\beta}\Vert\right\rbrace$$In some cases, instead of global optimization, it is sufficient to consider optimization by coordinates (widely studied in Daubechies et al. (2004)). If $f:\mathbb{R}^d\rightarrow\mathbf{R}$ is convex and differentiable, if $\mathbf{x}$ satisfies $f(\mathbf{x}+h\boldsymbol{e}_i)\geq f(\mathbf{x})$ for any $h>0$ and $i\in\{1,\cdots, d\}$then $f(\mathbf{x})=\min\{f\}$, where $\mathbf{e}=(\mathbf{e}_i)$ is the canonical basis of $\mathbb{R}^d$. However, this property is not true in the non-differentiable case. But if we assume that the non-differentiable part is separable (additively), it becomes true again. More specifically, if$$f(\mathbf{x})=g(\mathbf{x})+\sum_{i=1}^d h_i(x_i)$$with$$\left\lbrace\begin{array}{l}g: \mathbb{R}^d\rightarrow\mathbb{R}\text{ convex-differentiable}\\h_i: \mathbb{R}\rightarrow\mathbb{R}\text{ convex}\end{array}\right.$$This was the case for Lasso regression, $\beta)\mapsto\| \mathbf{y}-\beta_0-\mathbf{X}\beta\|_{\ell_2 }+\lambda\|\beta\|_{\ell_1}$, as shown by Tsen (2001). Getting back to our initial notations, we can use a coordinate descent algorithm: from an initial value $\mathbf{x}^{(0)}$, we consider (by iterating)$$x_j^{(k)}\in\text{argmin}\big\lbrace f(x_1^{(k)},\cdots,x_{k-1}^{(k)},x_k,x_{k+1}^{(k-1)},\cdots,x_n^{(k-1)})\big\rbrace$$ for $j=1,2,\cdots,n$These algorithmic problems and numerical issues may seem secondary to econometricians. However, they are essential in automatic learning: a technique is interesting if there is a stable and fast algorithm, which allows to obtain a solution. These optimization techniques can be transposed: for example, this coordinate descent technique can be used in the case of SVM methods (known as “vector support” methods) when the space is not linearly separable, and the classification error must be penalized (we will come back to this technique in the next section).

## In-sample, out-of-sample and cross-validation

These techniques seem intellectually interesting, but we have not yet discussed the choice of the penalty parameter $\lambda$. But this problem is actually more general, because comparing two parameters $\widehat{\beta}_{\lambda_1}$ and $\widehat{\beta}_{\lambda_2}$ is actually comparing two models. In particular, if we use a Lasso method, with different thresholds $\lambda$, we compare models that do not have the same dimension. Previously, we have addressed the problem of model comparison from an econometric perspective (by penalizing overly complex models). In the learning literature, judging the quality of a model on the data used to construct it does not make it possible to know how the model will behave on new data. This is the so-called “generalization” problem. The traditional approach then consists in separating the sample (size $n$) into two parts: a part that will be used to train the model (the training database, in-sample, size $m$) and a part that will be used to test the model (the testing database, out-of-sample, size $n-m$). The latter then makes it possible to measure a real predictive risk. Suppose that the data are generated by a linear model $y_i=\mathbf{x}_i^T \beta_0+\varepsilon_i$ where $\varepsilon_i$ are independent and centred law achievements. The empirical quadratic risk in-sample is here$$\frac{1}{m}\sum_{i=1}^m\mathbb{E}\big([\mathbf{x}_i^T \widehat{\beta}-\mathbf{x}_i^T \beta_0]^2\big)=\mathbb{E}\big([\mathbf{x}_i^T \widehat{\beta}-\mathbf{x}_i^T \beta_0]^2\big),$$for any observation $i$. Assuming the residuals $\varepsilon$ Gaussian, then we can show that this risk is worth $\sigma^2 \text{trace} (\Pi_X)/m$ is $\sigma^2 p/m$. On the other hand, the empirical out-of-sample quadratic risk is here $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big)$$where $\mathbf{x}$ is a new observation, independent of the others. It can be noted that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big\vert \mathbf{x}\big)=\text{Var}\big(\mathbf{x}^T \widehat{\beta}\big\vert \mathbf{x}\big)=\sigma^2\mathbf{x}^T(\mathbf{x}^T\mathbf{x})^{-1}\mathbf{x},$$and by integrating with respect to $\mathbf{x}$, $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T\beta_0]^2\big)=\sigma^2\text{trace}\big(\mathbb{E}[\mathbf{x}\mathbf{x}^T]\mathbb{E}\big[(\mathbf{x}^T\mathbf{x})^{-1}\big]\big).$$The expression is then different from that obtained in-sample, and using the Groves & Rothenberg (1969) increase, we can show that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big) \geq \sigma^2\frac{p}{m}$$which is pretty intuitive, when we start thinking about it. Except in some simple cases, there is no simple (explicit) formula. Note, however, that if $\mathbf{X}\sim\mathcal{N}(0,\sigma^2 \mathbb{I})$, then $\mathbf{x}^T \mathbf{x}$ follows a Wishart law, and it can be shown that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big)=\sigma^2\frac{p}{m-p-1}.$$If we now look at the empirical version: if $\widehat{\beta}$ is estimated on the first $m$ observations,$$\widehat{\mathcal{R}}^{~\text{ IS}}=\sum_{i=1}^m [y_i-\boldsymbol{x}_i^T\widehat{\boldsymbol{\beta}}]^2\text{ and }\widehat{\mathcal{R}}^{\text{ OS}}=\sum_{i=m+1}^{n} [y_i-\boldsymbol{x}_i^T\widehat{\boldsymbol{\beta}}]^2$$and as Leeb (2008) noted, $\widehat{\mathcal{R}}^{\text{IS}}-\widehat{\mathcal{R}}^{\text{OS}}\approx 2\cdot\nu$ where $\nu$ represents the number of degrees of freedom, which is not unlike the penalty used in the Akaike test.

Figure 4 shows the respective evolution of $\widehat{\mathcal{R}}^{\text{IS}}$ and $\widehat{\mathcal{R}}^{\text{OS}}$ according to the complexity of the model (number of degrees in a polynomial regression, number of nodes in splines, etc). The more complex the model, the more $\widehat{\mathcal{R}}^{\text{IS}}$ will decrease (this is the red curve, below). But that’s not what we’re interested in here: we want a model that predicts well on new data (i. e. out-of-sample). As Figure 4 shows, if the model is too simple, it does not predict well (as it does with in-sample data). But what we can see is that if the model is too complex, we are in a situation of “overlearning”: the model will start to model the noise. Of course, this figure should remind us of the one we’ve seen in our second post of that series

Figure 4 : Generalization, under- and over-fitting

Instead of splitting the database in two, with some of the data that will be used to calibrate the model and some to study its performance, it is also possible to use cross-validation. To present the general idea, we can go back to the “jackknife”, introduced by Quenouille (1949) (and formalized by Quenouille (1956) and Tukey (1958)) relatively used in statistics to reduce bias. Indeed, if we assume that $\{y_1,\cdots,y_n\}$ is a sample drawn according to a law $F_\theta$, and that we have an estimator $T_n (\mathbf{y})=T_n (y_1,\cdots,y_n)$, but that this estimator is biased, with $\mathbf{E}[T_n (\mathbf{Y})]=\theta+O(n^{-1})$, it is possible to reduce the bias by considering $$\widetilde{T}_n(\mathbf{y})=\frac{1}{n}\sum_{i=1}^n T_{n-1}(\mathbf{y}_{(i)})\text{ where }\mathbf{y}_{(i)}=(y_1,\cdots,y_{i-1},y_{i+1},\cdots,y_n)$$It can then be shown that $\mathbb{E}[\tilde{T}_n(Y)]=\theta+O(n^{-2})$The idea of cross-validation is based on the idea of building an estimator by removing an observation. Since we want to build a predictive model, we will compare the forecast obtained with the estimated model, and the missing observation$$\widehat{\mathcal{R}}^{\text{ CV}}=\frac{1}{n}\sum_{i=1}^n \ell(y_i,\widehat{m}_{(i)}(\mathbf{x}_i))$$We will speak here of the “leave-one-out” (loocv) method.

This technique reminds us of the traditional method used to find the optimal parameter in exponential smoothing methods for time series. In simple smoothing, we will construct a forecast from a time series as ${}_t\widehat{y}_{t+1} =\alpha\cdot{}_{t-1}\widehat{y}_t +(1-\alpha)\cdot y_t$, where $\alpha\in[0,1]$, and we will consider as “optimal” $$\alpha^\star = \underset{\alpha\in[0,1]}{\text{argmin}}\left\lbrace \sum_{t=2}^T \ell({}_{t-1}\widehat{y}_{t},y_{t}) \right\rbrace$$as described by Hyndman et al (2009).

The main problem with the leave-one-out method is that it requires calibration of n models, which can be problematic in large dimensions. An alternative method is cross validation by $k$-blocks (called “$k$-fold cross validation”) which consists in using a partition of $\{1,\cdots,n\}$ in $k$ groups (or blocks) of the same size, $\mathcal{I}_1,\cdots,\mathcal{I}_k$, and let us note $\mathcal{I}_{\bar j}=\{1,\cdots,n\}\setminus \mathcal{I}_j$. By noting $\widehat{m}_{(j)}$ built on the sample $\mathcal{I}_{\bar j}$, we then set:$$\widehat{\mathcal{R}}^{k-\text{ CV}}=\frac{1}{k}\sum_{j=1}^k \mathcal{R}_j\text{ where }\mathcal{R}_j=\frac{k}{n}\sum_{i\in\mathcal{I}_{{j}}} \ell(y_i,\widehat{m}_{(j)}(\mathbf{x}_i))$$Standard cross-validation, where only one observation is removed each time (loocv), is a special case, with $k=n$. Using $k=5$ or $10$ has a double advantage over $k=n$: (1) the number of estimates to be made is much smaller, 5 or 10 rather than $n$; (2) the samples used for estimation are less similar and therefore less correlated to each other, which tends to avoid excess variance, as recalled by James et al. (2013).

Another alternative is to use boosted samples. Let $\mathcal{I}_b$ be a sample of size $n$ obtained by drawing with replacement in $\{1,\cdots,n\}$ to know which observations $(y_i,\mathbf{x}_i)$ will be kept in the learning population (at each draw). Note $\mathcal{I}_{\bar b}=\{1,\cdots,n\}\setminus\mathcal{I}_b$. By noting $\widehat{m}_{(b)}$ built on sample $\mathcal{I}_b$, we then set :$$\widehat{\mathcal{R}}^{\text{ B}}=\frac{1}{B}\sum_{b=1}^B \mathcal{R}_b\text{ where }\mathcal{R}_b=\frac{n_{\overline{b}}}{n}\sum_{i\in\mathcal{I}_{\overline{b}}} \ell(y_i,\widehat{m}_{(b)}(\mathbf{x}_i))$$where $n_{\bar b}$ is the number of observations that have not been kept in $\mathcal{I}_b$. It should be noted that with this technique, on average $e^{-1}\sim36.7\%$ of the observations do not appear in the boosted sample, and we find an order of magnitude of the proportions used when creating a calibration sample, and a test sample. In fact, as Stone (1977) had shown, the minimization of AIC is to be compared to the cross-validation criterion, and Shao (1997) showed that the minimization of BIC corresponds to $k$-fold cross-validation, with $k=n/\log n$.

All those techniques here are mentioned in the “machine learning” section since they rely on automatic, computational techniques, and no probabilistic foundations are necessary. In many cases we did use the notation $m^\star$ (at least in the first posts on “machine learning” techniques) to highlight the fact that we want some sort of “optimal” model – and to make a distinction with estimators $\widehat{m}$ considered earlier, when we had some probabilistic framework. But of course, it is possible (and necessary) to build bridges between those two cultures…

References are online here. As explained in the introduction, it is some sort of online version of an introduction to our joint paper with Emmanuel Flachaire and Antoine Ly, Econometrics and Machine Learning (initially writen in French), that will actually appear soon in the journal Economics and Statistics (in English and in French).

# Foundations of Machine Learning, part 1

This post is the fifth one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 4 is online here.

In parallel with these tools developed by, and for economists, a whole literature has been developed on similar issues, centered on the problems of prediction and forecasting. For Breiman (2001a), a first difference comes from the fact that the statistic has developed around the principle of inference (or to explain the relationship linking $y$ to variables $\mathbf{x}$) while another culture is primarily interested in prediction. In a discussion that follows the article, David Cox states very clearly that in statistic (and econometrics) “predictive success (…) is not the primary basis for model choice“. We will get back here on the roots of automatic learning techniques. The important point, as we will see, is that the main concern of machine learning is related to the generalization properties of a model, i.e. its performance – according to a criterion chosen a priori – on new data, and therefore on non-sample tests.

## A learning machine

Today, we speak of “machine learning” to describe a whole set of techniques, often computational, as alternatives to the classical econometric approach. Before characterizing them as much as possible, it should be noted that historically other names have been given. For example, Friedman (1997) proposes to make the link between statistics (which closely resemble econometric techniques – hypothesis testing, ANOVA, linear regression, logistics, GLM, etc.) and what was then called “data mining” (which then included decision trees, methods from the closest neighbours, neural networks, etc.). The bridge between those two cultures corresponds to “statistical learning” techniques described in Hastie et al (2009). But one should keep in mind that machine learning is a very large field of research.

The so-called “natural” learning (as opposed to machine learning) is that of children, who learn to speak, read and play. Learning to speak means segmenting and categorizing sounds, and associating them with meanings. A child also learns simultaneously the structure of his or her mother tongue and acquires a set of words describing the world around him or her. Several techniques are possible, ranging from rote learning, generalization, discovery, more or less supervised or autonomous learning, etc. The idea in artificial intelligence is to take inspiration from the functioning of the brain to learn, to allow “artificial” or “automatic” learning, by a machine. A first application was to teach a machine to play a game (tic-tac-toe, chess, go, etc.). An essential step is to explain the objective it must achieve to win. One historical approach has been to teach the machine the rules of the game. If it allows you to play, it will not help the machine to play well. Assuming that the machine knows the rules of the game, and that it has a choice between several dozen possible moves, which one should it choose? The classical approach in artificial intelligence uses the so-called min-max algorithm using an evaluation function: in this algorithm, the machine searches forward in the possible moves tree, as far as the calculation resources allow (about ten moves in chess, for example). Then, it calculates different criteria (which have been previously indicated to her) for all positions (number of pieces taken, or lost, occupancy of the center, etc. in our example of the chess game), and finally, the machine plays the move that allows it to maximize its gain. Another example may be the classification and recognition of images or shapes. For example, the machine must identify a number in a handwritten handwriting (checks, ZIP code on envelopes, etc). It is a question of predicting the value of a variable $y$, knowing that a priori $y\in\{0,1,2,\cdots,8,9\}$. A classical strategy is to provide the machine with learning bases, in other words here millions of labelled (identified) images of handwritten numbers. A simple (and natural) strategy is to use a decision criterion based on the closest neighbors whose labels are known (using a predefined metric).

The method of the closest neighbors (“$k$-nearest neighbors”) can be described as follows: we consider (as in the previous part) a set of n observations, i. e. pairs $(y_i,\mathbf{x}_i)$ with $\mathbf{x}_i\in\mathbb{R}^p$. Let us consider a distance $\Delta$ on $\mathbb{R}^p$ (the Euclidean distance or the Mahalanobis distance, for example). Given a new observation $\mathbf{x}\in\mathbb{R}^p$, let us assume the ordered observations as a function of the distance between the $\mathbf{x}_i$ and $\mathbf{x}$, in the sense that $$\Delta(\mathbf{x}_1, \mathbf{x})\leq\Delta(\mathbf{x}_2, \mathbf{x})\leq\cdots\leq\Delta(\mathbf{x}_n, \mathbf{x})$$ then we can consider as prediction for y the average of the nearest $k$ neighbours,$$\widehat{m}_k(\mathbf{x})=\frac{1}{k}\sum_{i=1}^k y_i$$Learning here works by induction, based on a sample (called the learning – or training – sample).

Automatic learning includes those algorithms that give computers the ability to learn without being explicitly programmed (as Arthur Samuel defined it in 1959). The machine will then explore the data with a specific objective (such as searching for the nearest neighbours in the example just described). Tom Mitchell proposed a more precise definition in 1998: a computer program is said to learn from experience $E$ in relation to a task $T$ and a performance measure $P$, if its performance on $T$, measured by $P$, improves with experience $E$. Task $T$ can be a defect score for example, and performance $P$ can be the percentage of errors made. The system learns if the percentage of predicted defects increases with experience.

As we can see, machine learning is basically a problem of optimizing a criterion based on data (from now on called learning). Many textbooks on machine learning techniques propose algorithms, without ever mentioning any probabilistic model. In Watt et al (2016) for example, the word “probability” is mentioned only once, with this footnote that will surprise and make smile any econometricians, “the logistic regression can also be interpreted from a probabilistic perspective” (page 86). But many recent books offer a review of machine learning approaches using probabilistic theories, following the work of Vaillant and Vapnik. By proposing the paradigm of “probably almost correct” learning (PAC), a probabilistic flavor has been added to the previously very computational approach, by quantifying the error of the learning algorithm (usually in a classification problem).

To be continued (references are online here)…

# L’IA pour prédire les émeutes ?

Il y a quelques semaines, j’avais été contacte par un journaliste qui voulait me poser des questions, suite a notre article Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media. Ça a été l’occasion de me replonger dedans… et de voir ce qui a été écrit depuis… Et ce soir, je découvre un peu par hasard que l’article est paru, dans le numéro de Février de Science & Vie…

# Probabilistic Foundations of Econometrics, part 3

This post is the third one of our series on the history and foundations of econometric and machine learning models. Part 2 is online here.

## Exponential family and linear models

The Gaussian linear model is a special case of a large family of linear models, obtained when the conditional distribution of $Y$ (given the covariates) belongs to the exponential family$$f(y_i|\theta_i,\phi)=\exp\left(\frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+c(y_i,\phi)\right)$$ with $\theta_i=\psi(\mathbf{x}_i^T \beta)$. Functions $a$, $b$ and $c$ are specified according to the type of exponential law (studied extensively in statistics since Darmoix (1935), as Brown (1986) reminds us), and $\psi$ is a one-to-one mapping that the user must specify. Log-likelihood then has a simple expression $$\log\mathcal{L}(\mathbf{\theta},\phi|\mathbf{y}) =\frac{\sum_{i=1}^ny_i\theta_i-\sum_{i=1}^nb(\theta_i)}{a(\phi)}+\sum_{i=1}^n c(y_i,\phi)$$ and the first order condition is then written$$\frac{\partial \log \mathcal{L}(\mathbf{\theta},\phi|\mathbf{y})}{\partial \mathbf{\beta}} = \mathbf{X}^T\mathbf{W}^{-1}[\mathbf{y}-\widehat{\mathbf{y}}]=\mathbf{0}$$based on Müller’s (2011) notations, where $\mathbf{W}$ is a weight matrix (which depends on $\beta$). Given the link between $\theta$ and the expectation of $Y$, instead of specifying the function $\psi(\cdot)$, we will tend to specify the link function $g(\cdot)$ defined by $$\widehat{y}=m(\mathbf{x})=\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=g^{-1} (\mathbf{x}^T \beta)$$For the Gaussian linear regression we consider an identity link, while for the Poisson regression, the natural link (called canonical) is the logarithmic link. Here, as $\mathbf{W}$ depends on $\beta$ (with $\mathbf{W}=diag(\nabla g(\widehat{\mathbf{y}})Var[\mathbf{y}])$ there is generally no explicit formula for the maximum likelihood estimator. But an iterative algorithm makes it possible to obtain a numerical approximation. By setting $$\mathbf{z}=g(\widehat{\mathbf{y}})+(\mathbf{y}-\widehat{\mathbf{y}})\cdot\nabla g(\widehat{\mathbf{y}})$$ corresponding to the error term of a Taylor development in order 1 of $g$, we obtain an algorithm of the form$$\widehat{\beta}_{k+1}=[\mathbf{X}^T \mathbf{W}_k^{-1} \mathbf{X}]^{-1} \mathbf{X}^T \mathbf{W}_k^{-1} \mathbf{z}_k$$By iterating, we will define $\widehat{\beta}=\widehat{\beta}_{\infty}$, and we can show that – with some additional technical assumptions (detailed in Müller (2011)) – this estimator is asymptotically Gaussian, with $$\sqrt{n}(\widehat{\beta} -\beta)\overset{\mathcal{L}}{\rightarrow} \mathcal{N}(\mathbf{0},I(β)^{-1})$$where numerically $I(\beta)=\varphi\cdot[\mathbf{X}^T \mathbf{W}_\infty^{-1} \mathbf{X}]$.

From a numerical point of view, the computer will solve the first-order condition, and actually, the law of $Y$ does not really intervene. For example, one can estimate a “Poisson regression” even when observations are not integers (but they need to be positive). In other words, the law of $Y$ is only an interpretation here, and the algorithm could be introduced in a different way (as we will see later on), without necessarily having an underlying probabilistic model.

## Logistic Regression

Logistic regression is the generalized linear model obtained with a Bernoulli’s law, and a link function which is the quantile function of a logistic law (which corresponds to the canonical link in the sense of the exponential family). Taking into account the form of Bernoulli’s law, econometrics proposes a model for $y_i\in\{0,1\}$, in which the logarithm of the odds follows a linear model: $$\log\left(\frac{\mathbb{P}[Y=1\vert \mathbf{X}=\mathbf{x}]}{\mathbb{P}[Y\neq 1\vert \mathbf{X}=\mathbf{x}]}\right)=\beta_0+\mathbf{x}^T\beta$$or $$\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=\frac{e^{\beta_0+\mathbf{x}^T\beta}}{1+ e^{\beta_0+\mathbf{x}^T\beta}}=H(\beta_0+\mathbf{x}^T\beta)$$where $H(\cdot)=\exp(\cdot)/(1+exp(\cdot))$ is the cumulative distribution function of the logistic law. The estimation of $(\beta_0,\beta)$ is performed by maximizing the likelihood: $$\mathcal{L}=\prod_{i=1}^n \left(\frac{e^{\mathbf{x}_i^T\mathbf{\beta}}}{1+e^{\boldsymbol{x}_i^T\mathbf{\beta}}}\right)^{y_i}\left(\frac{1}{1+e^{\mathbf{x}_i^T\mathbf{\beta}}}\right)^{1-y_i}$$ It is said to be a linear models because isoprobability curves here are the parallel hyperplanes $b+\mathbf{x}^T\beta$. Rather than this model, popularized by Berkson (1944), some will prefer the probit model (see Berkson, 1951), introduced by Bliss (1934). In this model: $$\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=\Phi (\beta_0+\mathbf{x}^T\beta)$$

where $\Phi$ denotes the distribution function of the reduced centred normal distribution. This model has the advantage of having a direct link with the Gaussian linear model, since $y_i=\mathbf{1}(y_i^\star>0)$ with $y_i^\star=\beta_0+\mathbf{x}^T \beta+\varepsilon_i$ where the residuals are Gaussian, $\mathcal{N}(0,\sigma^2)$. An alternative is to have centered residuals of unit variance, and to consider a latent modeling of the form $y_i=\mathbf{1}(y_i^\star>\xi)$ (where $\xi$ will be fixed). As we can see, these techniques are fundamentally linked to an underlying stochastic model. In the body of the article, we present several alternative techniques – from the learning literature – for this classification problem (with two classes, here $0$ and $1$).

## Regression in high dimension

As we mentioned earlier, the first order condition $\mathbf{X}^T (\mathbf{X}\widehat{\beta}-\mathbf{y})=\mathbf{0}$ is solved numerically by performing a QR decomposition, at a cost which consists in $O(np^2)$ operations (where $p$ is the rank of $\mathbf{X}^T \mathbf{X}$). Numerically, this calculation can be long (either because $p$ is large or because $n$ is large), and a simpler strategy may be to sub-sample. Let $n_s\ll n$, and consider a sub-sample size $n_s$ of $\{1,\cdots,n\}$. Then $\widehat{\beta}_s=(\mathbf{X}_s^T \mathbf{X}_s )^{-1} \mathbf{X}_s^T\mathbf{y}_s$ is a good approximation of $\beta$ as shown by Dhillon et al. (2014). However, this algorithm is dangerous if some points have a high leverage (i.e. $L_i=\mathbf{x}_i(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i^T$). Tropp (2011) proposes to transform the data (in a linear way), but a more popular approach is to do non-uniform sub-sampling, with a probability related to the influence of observations (defined by $I_i=\widehat{\varepsilon}_iL_i/(1-L_i)^2$, and which unfortunately can only be calculated once the model is estimated).

In general, we will talk about massive data when the data table of size does not fit in the RAM memory of the computer. This situation is often encountered in statistical learning nowadays with very often $p\ll n$. This is why, in practice, many libraries of algorithms assimilated to machine learning use iterative methods to solve the first-order condition. When the parametric model to be calibrated is indeed convex and semi-differentiable, it is possible to use, for example, the stochastic gradient descent method as suggested by Bottou (2010). This last one allows to free oneself at each iteration from the calculation of the gradient on each observation of our learning base. Rather than making an average descent at each iteration, we start by drawing (without replacement) an observation $\mathbf{x}_i$ among the $n$ available. The model parameters are then corrected so that the prediction made from $\mathbf{x}_i$ is as close as possible to the true value $y_i$. The method is then repeated until all the data have been reviewed. In this algorithm there is therefore as much iteration as there are observations. Unlike the gradient descent algorithm (or Newton’s method) at each iteration, only one gradient vector is calculated (and no longer $n$). However, it is sometimes necessary to run this algorithm several times to increase the convergence of the model parameters. If the objective is, for example, to minimize a loss function $\ell$ between the estimator $m_\beta (\mathbf{x})$ and $y$ (like the quadratic loss function, as in the Gaussian linear regression) the algorithm can be summarized as follows:

• Step 0: Mix the data
• Iteration step: For $t=1,\cdots, n$, we pull $i\in\{1,\cdots,n\}$ without replacement, and we set $$\beta^{t+1} = \beta^{t} - \gamma_t\frac{ \partial{\ell(y_i,m_{\beta^t}(X_i)) } }{ \partial{ \beta}}$$

This algorithm can be repeated several times as a whole depending on the user’s needs. The advantage of this method is that at each iteration, it is not necessary to calculate the gradient on all observations (more sum). It is therefore suitable for large databases. This algorithm is based on a convergence in probability towards a neighborhood of the optimum (and not the optimum itself).

(references will be given in the very last post of that series) To be continued

# Probabilistic Foundations of Econometrics, part 2

This post is the second one of our series on the history and foundations of econometric and machine learning models. Part 1 is online here.

## Geometric Properties of this Linear Model

Let’s define the scalar product in $\mathbb{R}^n$, $⟨\mathbf{a},\mathbf{b}⟩=\mathbf{a}^T\mathbf{b}$, and let’s note $\|\cdot\|$ the associated Euclidean standard, $\|\mathbf{a}\|=\sqrt{\mathbf{a}^T\mathbf{a}}$ (denoted $\|\cdot\|_{\ell_2}$ in the next post). Note $\mathcal{E}_X$ the space generated by all linear combinations of the $\mathbf{X}$ components (adding the constant). If the explanatory variables are linearly independent, $\mathbf{X}$ is a full (column) rank matrix and $\mathcal{E}_X$ is a space of dimension $p+1$. Let’s assume from now on that the variables $\mathbf{x}$  and $y$ are centered here. Note that no law hypothesis is made in this section, the geometric properties are derived from the properties of expectation and variance in the set of finite variance variables.

With this notation, it should be noted that the linear model is written $m(\mathbf{x})=⟨\mathbf{x},\beta⟩$. The space $H_z=\{\mathbf{x}\in\mathbb{R}^{p+1}:m(\mathbf{x})=z\}$ is a hyperplane (affine) that separates the space in two. Let’s define the orthogonal projection operator on $\mathcal{E}_X$, $\Pi_X =\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T$. Thus, the forecast that can be made for it is: $$\widehat{\mathbf{y}}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T\mathbf{y}=\Pi_X\mathbf{y}$$. As, $\widehat{\varepsilon}=\mathbf{y}-\widehat{\mathbf{y}}=(\mathbb{I}-\Pi_X)\mathbf{y}=\Pi_{X^\perp}\mathbf{y}$, we note that $\widehat{\varepsilon}\perp\mathbf{x}$, which will be interpreted as meaning that residuals are a term of innovation, unpredictable in the sense that $\Pi_{X }\widehat{\varepsilon}=\mathbf{0}$. The Pythagorean theorem is written here: $$\Vert \mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y} \Vert^2+\Vert \Pi_{ {X}^\perp}\mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y}\Vert^2+\Vert \mathbf{y}-\Pi_{ {X}}\mathbf{y}\Vert^2=\Vert\widehat{\mathbf{y}}\Vert^2+\Vert\widehat{\mathbf{\varepsilon}}\Vert^2$$which is classically translated in terms of the sum of squares: $$\underbrace{\sum_{i=1}^n y_i^2}_{n\times\text{total variance}}=\underbrace{\sum_{i=1}^n \widehat{y}_i^2}_{n\times\text{explained variance}}+\underbrace{\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{n\times\text{residual variance}}$$The coefficient of determination, $R^2$, is then interpreted as the square of the cosine of the angle $\theta$ between $\mathbf{y}$ and $\Pi_X \mathbf{y}$ : $$R^2=\frac{\Vert \Pi_{{X}} \mathbf{y}\Vert^2}{\Vert \mathbf{y}\Vert^2}=1-\frac{\Vert \Pi_{ {X}^\perp} \mathbf{y}\Vert^2}{\Vert \mathbf {y}\Vert^2}=\cos^2(\theta)$$An important application was obtained by Frish & Waugh (1933), when the explanatory variables are divided into two groups, $\mathbf{X}=[\mathbf{X}_1 |\mathbf{X}_2]$, so that the regression becomes $y=\beta_0+\mathbf{X}_1 β_1+\mathbf{X}_2 β_2+\varepsilon$. Frish & Waugh (1933) showed that two successive projections could be considered. Indeed, if $\mathbf{y}_2^\star=\Pi_{X_1^\perp} \mathbf{y}$ and $X_2^\star=\Pi_{X_1^\perp}\mathbf{X}_2$, we can show that $$\widehat{\beta} _2=[{\mathbf{X}_2^\star}^T \mathbf{X}_2^\star]^{-1}{\mathbf{X}_2^\star}^T \mathbf{y}_2^\star$$ In other words, the overall estimate is equivalent to the combination of independent estimates of the two models if $\mathbf{X}_2^\star=\mathbf{X}_2$, i.e. $\mathbf{X}_2\in \mathcal{E}_{X_1}^\perp$, which can be noted $\mathbf{x}_1\perp\mathbf{x}_2$ We obtain here the Frisch-Waugh theorem which guarantees that if the explanatory variables between the two groups are orthogonal, then the overall estimate is equivalent to two independent regressions, on each of the sets of explanatory variables. This is a theorem of double projection, on orthogonal spaces. Many results and interpretations are obtained through geometric interpretations (fundamentally related to the links between conditional expectation and the orthogonal projection in space of variables of finite variance).

This geometric interpretation might help to get a better understanding of the problem of under-identification, i.e. the case where the real model would be $y_i=\beta_0+ \mathbf{x}_1^T \beta_1+\mathbf{x}_2^T \beta_2+\varepsilon_i$, but the estimated model is $y_i=b_0+\mathbf{x}_1^T \mathbf{b}_1+\eta_i$. The maximum likelihood estimator of $\mathbf{b}_1$ is $$\widehat{\mathbf{b}}_1=\mathbf {\beta}_1 + \underbrace{ (\mathbf {X}_1^T\mathbf {X}_1)^{-1} \mathbf {X}_1^T \mathbf {X}_{2} \mathbf{\beta}_2}_{\mathbf{\beta}_{12}}+\underbrace{(\mathbf{X}_1^{T}\mathbf{X}_1)^{-1} \mathbf{X}_1^T\varepsilon}_{\nu}$$so that $\mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1+\beta_{12}$, the bias ($\beta_{12}$) being null only in the case where $\mathbf{X}_1^T \mathbf{X}_2=\mathbf{0}$ (i. e. $\mathbf{X}_1\perp \mathbf{X}_2$): we find here a consequence of the Frisch-Waugh theorem.

On the other hand, over-identification corresponds to the case where the real model would be $y_i=\beta_0+\mathbf{x}_1^T \beta_1+\varepsilon_i$, but the estimated model is $y_i=b_0+ \mathbf{x}_1^T \mathbf{b} _1+\mathbf{x}_2^T \mathbf{b}_2+\eta_i$. In this case, the estimate is unbiased, in the sense that $\mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1$ but the estimator is not efficient. Later on, we will discuss an effective method for selecting variables (and avoid over-identification).

## From parametric to non-parametric

We can rewrite equation (4) in the form $\widehat{\mathbf{y}}=\Pi_X\mathbf{y}$ which helps us to see the forecast directly as a linear transformation of the observations. More generally, a linear predictor can be obtained by considering $m(\mathbf{x})=\mathbf{s}_{\mathbf{x}}^T \mathbf{y}$, where $\mathbf{s}_{\mathbf{x}}$ is a weight vector, which depends on $\mathbf{x}$, interpreted as a smoothing vector. Using the vectors $\mathbf{s}_{\mathbf{x}_i}$, calculated from the observations $\mathbf{x}_i$, we obtain a matrix $\mathbf{S}$ of size $n\times n$, and $\widehat{\mathbf{y}}=\mathbf{S}\mathbf{y}$. In the case of the linear regression described above, $\mathbf{s}_{\mathbf{x}}=\mathbf{X}[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{x}$, and in that case $\text{trace}(\mathbf{S})$ is the number of columns in the $\mathbf{X}$ matrix (the number of explanatory variables). In this context of more general linear predictors, $\text{trace}(\mathbf{S})$ is often seen as equivalent to the number of parameters (or complexity, or dimension, of the model), and $\nu=n-\text{trace}(\mathbf{S})$ is then the number of degrees of freedom (see Ruppert et al., 2003; Simonoff, 1996). The principle of parsimony says that we should minimize this dimension (the trace of the matrix $\mathbf{S}$) as much as possible. But in the general case, this dimension is more to obtain, explicitely.

The estimator introduced by Nadaraya (1964) and Watson (1964), in the case of a simple non-parametric regression, is also written in this form since$$\widehat{m}_h(x)=\mathbf{s}_{x}^T\mathbf{y}=\sum_{i=1}^n \mathbf{s}_{x,i}y_i$$where$$\mathbf{s}_{x,i}=\frac{K_h(x-x_i)}{K_h(x-x_1)+\cdots+K_h(x-x_n)}$$ where $K(\cdot)$ is a kernel function, which assigns a value that is lower the closer $x_i$ is to $x$, and $h>0$ is the bandwidth. The introduction of this metaparameter $h$ is an important issue, as it should be chosen wisely. Using asymptotic developments, we can show that if $X$ has density $f$, $$\text{biais}[\widehat{m}_h(x)]=\mathbb{E}[\widehat{m}_h(x)]-m(x)\sim {h^2}\left(\frac{C_1 }{2}m''(x)+C_2 m'(x)\frac{f'(x)}{f(x)}\right)$$and $$\displaystyle{{\text{Var}[\widehat{m}_h(x)]\sim\frac{C_3}{{nh}}\frac{\sigma(x)}{f(x)}}}$$for some constants that can be estimated (see Simonoff (1996) for a discussion). These two functions evolve inversely with $h$, as shown in Figure 1 (where the metaparameter on the $x$-axis is here, actually, $h^{-1}$). Keep in ming that we will see a similar graph in the context of machine learning models.

Figure 1. Choice of meta-parameter and the Goldilocks problem: it must not be too large (otherwise there is too much variance), nor too small (otherwise there is too much bias).

The natural idea is then to try to minimize the mean square error, the MSE, defined as $bias[\widehat{m}_h (x)]^2+Var[\widehat{m}_h (x)]$, and them integrate over $x$, which gives an optimal value for $h$ of the form $h^\star=O(n^{-1/5})$, and reminds us of Silverman’s rule – see Silverman (1986). In larger dimensions, for continuous $\mathbf{x}$ variables, a multivariate kernel with matrix bandwidth $\mathbf{H}$ can be used, and$$\mathbb{E}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim m(\mathbf{x})+\frac{C_1}{2}\text{trace}\big(\mathbf{H}^Tm''(\mathbf{x})\mathbf{H}\big)+C_2\frac{m'(\boldsymbol{x})^T\mathbf{H}\mathbf{H}^T \nabla f(\mathbf{x})}{f(\mathbf{x})}$$while$$\text{Var}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim\frac{C_3}{n~\text{det}(\mathbf{H})}\frac{\sigma(\mathbf{x})}{f(\mathbf{x})}$$
If $\mathbf{H}$ is a diagonal matrix, with the same term $h$  on the diagonal, then $h^\star=O(n^{-1/(4+dim(\mathbf{x}))}$. However, in practice, there will be more interest in the integrated version of the quadratic error, $$MISE(\widehat{m}_{h})=\mathbb{E}[MSE(\widehat{m}_{h}(X))]=\int MSE(\widehat{m}_{h}(x))dF(x)$$and we can prove that $$MISE[\widehat{m}_h]\sim \overbrace{\frac{h^4}{4}\left(\int x^2k(x)dx\right)^2\int\big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]^2dx}^{\text{bias}^2} +\overbrace{\frac{\sigma^2}{nh}\int k^2(x)dx \cdot\int\frac{dx}{f(x)}}^{\text{variance}}$$as n→∞ and nh→∞. Here we find an asymptotic relationship that again recalls Silverman’s (1986) order of magnitude, $$h^\star =n^{-\frac{1}{5}}\left(\frac{C_1\int \frac{dx}{f(x)}}{C_2\int \big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]dx}\right)^{\frac{1}{5}}$$The main problem here, in practice, is that many of the terms in the expression above are unknown. Automatic learning offers computational techniques, when the econometrician used to searching for asymptotic (mathematical) properties.

To be continued (references mentioned above are online here)…

# Probabilistic Foundations of Econometrics, part 1

In a series of posts, I wanted to get into details of the history and foundations of econometric and machine learning models. It will be some sort of online version of our joint paper with Emmanuel Flachaire and Antoine Ly, Econometrics and Machine Learning (initially writen in French), that will actually appear soon in the journal Economics and Statistics. This is the first one…

The importance of probabilistic models in economics is rooted in Working’s (1927) questions and the attempts to answer them in Tinbergen’s two volumes (1939). The latter have subsequently generated a great deal of work, as recalled by Duo (1993) in his book on the foundations of econometrics, and more particularly in the first chapter “The Probability Foundations of Econometrics”. It should be recalled that Trygve Haavelmo was awarded the Nobel Prize in Economics in 1989 for his “clarification of the foundations of the probabilistic theory of econometrics”. Because as Haavelmo (1944) (initiating a profound change in econometric theory in the 1930s, as recalled in Morgan’s Chapter 8 (1990)) showed, econometrics is fundamentally based on a probabilistic model, for two main reasons. First, the use of statistical quantities (or “measures”) such as means, standard errors and correlation coefficients for inferential purposes can only be justified if the process generating the data can be expressed in terms of a probabilistic model. Second, the probability approach is relatively general, and is particularly well suited to the analysis of “dependent” and “non-homogeneous” observations, as they are often found on economic data.We will then assume that there is a probabilistic space $(\Omega,\mathcal{F},\mathbb{P})$ such that observations $(y_i,\mathbf{x}_i)$ are seen as realizations of random variables $(Y_i, \mathbf{X}_i)$. In practice, however, we are not very interested in the joint law of the couple $(Y, \mathbf{X})$ : the law of $\mathbf{X}$ is unknown, and it is the law of Y conditional on $\mathbf{X}$ that will be interested in. In the following, we will note $x$ a single observation, $\mathbf{x}$ a vector of observations, $X$ a random variable, and $\mathbf{X}$ a random vector. Abusively, $\mathbf{X}$ may also designate the matrix of individual observations (denoted $\mathbf{x}_i$), depending on the context.

## Foundations of mathematical statistics

As recalled in Vapnik’s (1998) introduction, inference in parametric statistics is based on the following belief: the statistician knows the problem to be analyzed well, in particular, he knows the physical law that generates the stochastic properties of the data, and the function to be found is written via a finite number of parameters[1]. To find these parameters, the maximum likelihood method is used. The purpose of the theory is to justify this approach (by discovering and describing its favorable properties). We will see that in learning, philosophy is very different, since we do not have a priori reliable information on the statistical law underlying the problem, nor even on the function we would like to approach (we will then propose methods to construct an approximation from the data at our disposal, as in (1998)). A “golden age” of parametric inference, from 1930 to 1960, laid the foundations for mathematical statistics, which can be found in all statistical textbooks, including today. As Vapnik (1998) states, the classical parametric paradigm is based on the following three beliefs:

1. To find a functional relationship from the data, the statistician is able to define a set of functions, linear in their parameters, that contain a good approximation of the desired function. The number of parameters describing this set is small.
2. The statistical law underlying the stochastic component of most real-life problems is the normal law. This belief has been supported by reference to the central limit theorem, which stipulates that under large conditions the sum of a large number of random variables is approximated by the normal law.
3. The maximum likelihood method is a good tool for estimating parameters.

In this section we will come back to the construction of the econometric paradigm, directly inspired by that of classical inferential statistics.

## Conditional laws and likelihood

Linear econometrics has been constructed under the assumption of individual data, which amounts to assuming independent variables $(Y_i, \mathbf{X}_i)$ (if it is possible to imagine temporal observations – then we would have a process $(Y_t, \mathbf{X}_t)$ – but we will not discuss time series here). More precisely, we will assume that, conditionally to the explanatory variables $\mathbf{X}_i$, the variables $Y_i$ are independent. We will also assume that these conditional laws remain in the same parametric family, but that the parameter is a function of $\mathbf{x}$. In the Gaussian linear model it is assumed that: $$(Y\vert \mathbf{X}=\mathbf{x})\overset{\mathcal{L}}{\sim}\mathcal{N}(\mu(\mathbf{x}),\sigma^2)~~~~ (1)$$where $\mu(\mathbf{x})=\beta_0+\mathbf{x}^T\mathbf{\beta}$ and $\mathbf{\beta}\in\mathbb{R}^{p}$.

It is usually called a ‘linear’ model since $\mathbb{E}[Y\vert \mathbf{X}=\mathbf{x}]=\beta_0+\mathbf{x}^T\mathbf{\beta}$ is a linear combination of covariates[2]. It is said to be a homoscedastic model if $Var[Y|\mathbf{X}=\mathbf{x}]=\sigma^2$, where $\sigma^2$ is a positive constant. To estimate the parameters, the traditional approach is to use the Maximum Likelihood estimator, as initially suggested by Ronald Fisher. In the case of the Gaussian linear model, log-likelihood is written:  $$\log\mathcal{L}(\beta_0, \mathbf{\beta},\sigma^2\vert \mathbf{y},\mathbf{x}) = -\frac{n}{2}\log[2\pi\sigma^2] - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-\beta_0-\mathbf{x}_i^T\mathbf{\beta})^2$$Note that the term on the right, measuring a distance between the data and the model, will be interpreted as deviance in generalized linear models. Then we will set: $$(\widehat{\beta}_0,\widehat{\mathbf{\beta}},\widehat{\sigma}^2)=\text{argmax}\left\lbrace\log\mathcal{L}(\beta_0, \mathbf{\beta},\sigma^2\vert \mathbf{y},\mathbf{x})\right\rbrace$$The maximum likelihood estimator is obtained by minimizing the sum of the error squares (the so-called “least squares” estimator) that we will find in the “machine learning” approach.

The first order conditions allow to find the normal equations, whose matrix writing is $\mathbf{X}^T[\mathbf{y}-\mathbf{X}\mathbf{\beta}]=\mathbf{0}$, which can also be written $(\mathbf{X}^T \mathbf{X})\mathbf{\beta}=\mathbf{X}^T \mathbf{y}$. If $\mathbf{X}$ is a full (column) rank matrix, then we find the classical estimator:$$\widehat{\mathbf{\beta}}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}=\mathbf{\beta}+(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^{-1}\mathbf{\varepsilon}~~~(2)$$using residual-based writing (as often in econometrics), $y=\mathbf{x}^T\mathbf{\beta}+\varepsilon$. Gauss Markov’s theorem ensures that this estimator is the unbiased linear estimator with minimum variance. It can then be shown that $\widehat{\mathbf{\beta}}\sim\mathcal{N}(\mathbf{\beta},\sigma^2(\mathbf{X}^T\mathbf{X})^{-1})$, and in particular, if we simply need the first two moments : $$\mathbb{E}[\widehat{\mathbf{\beta}}]=\mathbf{\beta}~~~Var[\widehat{\mathbf{\beta}}]=\sigma^2 [\mathbf{X}^T\mathbf{X}]^{-1}$$In fact, the normality hypothesis makes it possible to make a link with mathematical statistics, but it is possible to construct this estimator given by equation (2) without that Gaussian assumption. Hence, if we assume that $Y|\mathbf{X}$ has the same distribution as $\mathbf{x}^T\mathbf{\beta}+\varepsilon$, where $\mathbb{E}[\varepsilon]=0$, $Var[\varepsilon]=\sigma^2$ and $Cov[X_j,\varepsilon]=0$ for all $j$, then $\widehat{\mathbf{\beta}}$ is an unbiased estimator of $\mathbf{\beta}$ with smallest variance[3] among unbiased linear estimators. Furthermore, if we cannot get normality at finite distance, asymptotically this estimator is Gaussian, with $$\sqrt{n}(\widehat{\mathbf{\beta}}-\mathbf{\beta})\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(\mathbf{0},\mathbf{\Sigma})$$as $n\rightarrow\infty$, for some matrix $\mathbf{\Sigma}$.
The condition of having a full rank $\mathbf{X}$ matrix can be (numerically) strong in large dimensions. If it is not satisfied, $(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T$ does not exist. If $\mathbb{I}$ denotes the identity matrix, however, it should be noted that $(\mathbf{X}^T \mathbf{X}+\lambda\mathbb{I})^{-1}\mathbf{X}^T$ still exists, whatever $\lambda>0$. This estimator is called the ridge estimator of level \lambda (introduced in the 1960s by Hoerl (1962), and associated with a regularization studied by Tikhonov (1963)). This estimator naturally appears in a Bayesian econometric context.

## Residuals

It is not uncommon to introduce the linear model from the distribution of the residuals, as we mentioned earlier. Also, equation (1) is written as often: $$y_i=\beta_0+\mathbf{x}_i^T\mathbf{\beta}+\varepsilon_i~~~~(3)$$where $\varepsilon_i$’s are realizations of independent and identically distributed random variables (i.i.d.) from some $\mathcal{N}(0,\sigma^2)$ distribution. With a vector notation, we will write $\mathbf{\varepsilon}\overset{\mathcal{L}}{\sim}\mathcal{N}(\mathbf{0},\sigma^2\mathbb{I})$. The estimated residuals are defined as: $$\widehat{\varepsilon}_i =y_i-[\widehat{\beta}_0+\mathbf{x}_i^T\widehat{\mathbf{\beta}}]$$ Those (estimated) residuals are basic tools for diagnosing the relevance of the model.

An extension of the model described by equation (1) has been proposed to take into account a possible heteroscedastic character: $$(Y\vert \mathbf{X}=\mathbf{x})\overset{\mathcal{L}}{\sim}\mathcal{N}(\mu(\mathbf{x}),\sigma^2(\mathbf{x}))$$where $\sigma^2(\mathbf{x})$ is a positive function of the explanatory variables. This model can be rewritten as: $$y_i=\beta_0+\mathbf{x}_i^T\mathbf{\beta}+\sigma^2(\mathbf{x}_i)\cdot\varepsilon_i$$where residuals are always i.i.d., with unit variance, $$\varepsilon_i=\frac{y_i-[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}{\sigma(\mathbf{x}_i)}$$ While residuals based equations are popular in linear econometrics (when the dependent variable is continuous), it is no longer popular in counting models, or logistic regression.

However, writing using an error term (as in equation (3)) raises many questions about the representation of an economic relationship between two quantities. For example, it can be assumed that there is a relationship (linear to begin with) between the quantities of a traded good, $q$ and its price $p$. This allows us to imagine a supply equation$$q_i=\beta_0+\beta_1 p_i+u_i$$($u_i$ being an error term) where the quantity sold depends on the price, but in an equally legitimate way, one can imagine that the price depends on the quantity produced (what one could call a demand equation), $$p_i=\alpha_0+\alpha_1 q_i+v_i$$($v_i$ denoting another error term). Historically, the error term in equation (3) could be interpreted as an idiosyncratic error on the variable $y$, the so-called explanatory variables being assumed to be fixed, but this interpretation often makes the link between an economic relationship and a complicated economic model difficult, the economic theory speaking abstractly about a relationship between a magnitude, the econometric model imposing a specific shape (what magnitude is $y$ and what magnitude is $x$) as shown in more detail in Morgan (1990) Chapter 7.

(references mentioned above are online here). To be continued…

[1] This approach can be compared to structural econometrics, as presented for example in Kean (2010).

[2] Here, we will try to distinguish $\beta_0$, the intercept, and the other parameters $\mathbf{\beta}$, since they are considered differently in many extensions (e.g. regularization). Nevertheless, in many expressions $\mathbf{\beta}$ will denote the joint vector $(\beta_0, \mathbf{\beta})$, for general formulas, to avoid too heavy notations.

[3] In the sense that the difference between variance matrices is a positive matrix.

# On my way to Québec City

This Friday, I will be in Quebec city to give a talk at the Actuarial Departement, on risk pooling and price segmentation. Slides are now available

# Gini Regressions and Heteroskedasticity

Our joint paper “Gini Regressions and Heteroskedasticity” with Ndéné Ka, Stéphane Mussard and Oumar Hamady Ndiaye just appear in Econometrics.

We propose an Aitken estimator for Gini regression. The suggested A-Gini estimator is proven to be a U-statistics. Monte Carlo simulations are provided to deal with heteroskedasticity and to make some comparisons between the generalized least squares and the Gini regression. A Gini-White test is proposed and shows that a better power is obtained compared with the usual White test when outlying observations contaminate the data

# Les réseaux pour réinventer l’assurance ?

La théorie des réseaux, ou des graphs, est née en 1735, suite aux travaux de Léonard Euler, qui essayait de trouver une promenade – à partir d’un point donné – qui fasse revenir à ce point en passant une fois et une seule par chacun des sept ponts de la ville de Königsberg. On peut rapprocher ces réseaux des réseaux de métro, constitués de stations (les nœuds), liés entre deux par des rails, ou pas, ou plus généralement un réseau routier, pouvant donner lieu à des études de congestion, par exemple. Mais les réseaux sont aujourd’hui surtout sociaux, reliant les personnes, par des liens d’amitiés, professionnels, familiaux, ou monétaires. L’analyse des réseaux permet de créer des communautés relativement homogène, acceptant de partager un risque, recréant une mutualisation.

# Réseau et crédit

En généalogie, on aura des réseaux hiérarchiques, un enfant étant lié à ses parents, eux-mêmes reliés à leurs parents. En sociologie, les réseaux sociaux permettent d’analyser les liens entre des individus (ou des organisations) au sein d’un ensemble. On pourra étudier les amitiés dans une cour d’école (un lien pouvant être une invitation à un anniversaire) ou des échanges de messages électroniques dans une entreprise (la base des courriels d’Enron a ainsi été abondamment utilisée, avec plus de 180 000 messages échangés entre 36 000 employés[i]). La Figure 1 montre ainsi deux réseaux de 20 individus (A, B, …, T).

Figure 1 : réseaux aléatoires, 20 nœuds (Watts-Strogatz et Barbasi),

Dans une vision de type Facebook ou Linkedin, on dira que E et F sont liés, au sens « amis », s’il existe un segment reliant les points E et F. Un réseau peut être dirigé, par exemple si on étudie les échanges de messages (E a écrit à F), ou des prêts d’argent (E a prêté de l’argent à F). Si historiquement seule l’adjacence était étudiée (existence ou non de liens), on peut aujourd’hui rajouter des poids, par exemple le montant d’un prêt financier. Babutsidze (2012) étudie ainsi les positions de banques françaises et allemandes dans les prêts interbancaires au sein de la zone Europe (les nœuds sont alors les banques). L’étude des réseaux au sein de communautés de villages dans les pays en développement a permis de mieux comprendre les mécanismes de finance informelle. Banerjee et al. (2013) étudient ainsi la diffusion de l’information dans un réseau, et plus particulièrement les prêts de microfinance.

Si les réseaux sont utiles pour mieux organiser le microcrédit, CNN notait en 2015 que Facebook permettait à des organismes de crédit d’utiliser le réseau social d’un emprunteur pour déterminer s’il représente un bon risque de crédit, ou pas. En particulier, si le score de crédit des amis était trop faible, une personne pouvait se voir refuser un crédit. Cette situation est dangereuse à cause de propriétés particulières des réseaux, et plus particulièrement le paradoxe des amis.

# Du tout petit monde au paradoxe des amis

En 1929, Frigyes Karinthy a émis l’hypothèse que toute personne sur terre pouvait être relliée à n’importe quelle autre par une succession de relations individuelles comprenant au plus 6 maillons. « Nous devrions sélectionner n’importe quelle personne du 1,5 milliard d’habitants de la planète, n’importe qui, n’importe où. Il paraît que, n’utilisant pas plus de cinq individus, l’un d’entre eux étant une connaissance personnelle, il pourrait contacter les individus choisis en n’utilisant rien d’autre que le réseau des connaissances personnelles ». Cette théorie des six poignées de main a vu son origine  dans une nouvelle littéraire. Il faudra attendre les travaux de

Michael Gurevich dans les années 60, puis Stanley Milgram dix ans après, pour voir apparaître les premières tentatives de quantification de ces relations, sous le nom de « Small World Problem ». Si Leskovec & Horvitz (2008) ont confirmé cet ordre de grandeur, via l’analyse de plusieurs milliards de messages échangés à l’aide de la plateforme Windows Live Messenger, plus récemment, Baghat et al. (2016) ont estimé que deux personnes quelconques sur Facebook étaient connectées par une moyenne de trois personnes et demi. Sur le réseau aléatoire de gauche, une personne a, en moyenne, 2 amis, alors qu’un ami pris au hasard a en moyenne 2.25 amis. Sur le réseau de droite, l’écart est encore plus important, car si là aussi une personne a, en moyenne, 2 amis, un ami pris au hasard aura en moyenne plus de 4 amis.

Figure 2 : réseaux aléatoires, 500 nœuds (Watts-Strogatz et Barbasi)

Ce paradoxe, observé en 1991 par le sociologue Scott Feld se démontre très facilement. Heuristiquement, on peut voir un lien avec la propriété probabiliste$$\frac{\mathbb{E}[X^2]}{\mathbb{E}[X]}=\mathbb{E}[X]+\frac{\text{Var}[X]}{\mathbb{E}[X]}>\mathbb{E}[X]$$où le terme de gauche est le nombre d’amis de mes amis, divisé par mon nombre d’amis. La différence est d’autant plus grande que la dispersion du nombre d’amis est importante. Si le réseau de gauche est très dense, celui de droite en revanche possède une propriété de loi puissance : la distribution du nombre d’amis suit une loi en fonction puissance (ou loi de Zipf, ou de Pareto). La Figure 3 montre la distribution du nombre d’amis sur un réseau, dans double échelle logarithmique : la linéarité indique une distribution en fonction puissance. On retrouve ce genre de distribution dans un très grand nombre de réseaux, en particulier Facebook, comme l’a montré Wohlgemuth & Matache (2014).

Figure 3 : distribution du nombre d’amis sur des réseaux aléatoires simulés (Watts-Strogatz et Barbasi en rouge)

L’interprétation classique est que certaines personnes sont centrales dans le réseau, avec un très grand nombre de connexions. Cette propriété est très connue en marketing (on parlera alors d’effet de pair, « peer effect ») mais elle a des impacts aussi en gestion des risques, ou en santé publique. Chrisakis & Fowler (2010) ont ainsi montré que les épidémies de grippe peuvent être détectées près de deux semaines en avance, en surveillant l’infection dans un réseau social. En particulier, l’analyse de la santé des personnes centrales dans un réseau est « an ideal way to predict outbreaks, but detailed information doesn’t exist for most groups, and to produce it would be time-consuming and costly ». Pour revenir à l’exemple du score de crédit, si ce dernier se trouve être corrélé aux nombres d’amis, le paradoxe des amis rend dangereuse l’utilisation du score des amis comme indicateur du risque d’un individu !

# L’importance de l’homophilie

Un autre trait important des réseaux est la notion d’homophilie, introduite en 2001 en sociologie par deux articles importants, correspondant à la tendance à être connecté à ses semblables. McPherson et al (2001) partait du principe que la similitude engendre la connexion, et par conséquent, les réseaux personnels des gens sont homogènes sur de nombreuses caractéristiques sociodémographiques, comportementales et intrapersonnelles. Moody (2001) étudiait ainsi les amitiés dans les cours de récréations à l’école primaire, aux États-Unis, et plus particulièrement les amitiés interraciales. Easley & Kleinberg (2010) présente ainsi de nombres conséquences de l’homophilie, allant de la constitution des tables lors de repas d’affaire, à l’attribution de crédit aux États-Unis. La mesure de l’homophilie revient à se demander, compte tenu de groupes préexistants (en fonction du genre, de l’âge, de la catégorie socioprofessionnelle, etc) comment se répartissent les liaisons, entre les groupes, ou à intérieur des groupes.

Figure 4 : faible homophilie (en haut) et forte homophilie (en bas).

Dans un contexte d’assurance, un actuaire cherche à créer des classes tarifaires, des groupes homogènes en termes de risques, suivant des variables explicatives (les variables dites tarifaires). Les personnes qui habitent au même endroit, qui conduisent les mêmes types de véhicules, et qui ont les mêmes caractéristiques, auront de fortes chances d’être dans la même classe. Mais si l’homophilie existe dans une population, un groupe tarifaire pourrait peut-être s’observer sur réseau d’amis. Pourquoi ne pas alors envisager de créer des groupes au sein d’un réseau ?

# Utiliser les réseaux en assurance

Dans cet esprit, en 2010, Friendsurance a été lancé en Allemagne et compte en 2018 plus de 100000 assurés[i]. En France, une courte expérience d’assurance collaborative avait été lancée en 2015, avec Inspeer, proposant de mutualiser avec ses proches les franchises d’assurance dommage (en assurance auto ou habitation) entre amis. Ces types d’assurances collaboratives, parfois appelés assurance pair à pair – ou « peer-to-peer insurance » –  reposent sur la constitution de petits groupes par un courtier. Une partie des primes d’assurance versées est versée à un fonds collectif, l’autre partie à une compagnie d’assurance tierce. Les dommages mineurs subis par le preneur d’assurance sont d’abord pris en charge par ce fonds de groupe. Pour les sinistres dépassant la franchise, il est fait appel à l’assureur habituel. Un groupe peut être constitué par les assurés, formant un réseau social un peu comme Facebook. Dans ce modèle, la seule exigence est que tous les membres du groupe doivent avoir le même type d’assurance (par exemple une assurance responsabilité civile avec une assurance des frais juridiques).

Comme le notait Schiller (2013), ce type de mécanisme a beaucoup de vertus, la première étant de diminuer les coûts, et le risque de fraude. On n’a en effet pas tendance à tricher sur le coût d’un sinistre lorsque le risque est porté par des membres de la famille, ou des amis. L’anonymat de la mutualité qui existe dans la loi des grands nombres disparait. Mais n’est-on pas en train de réinventer la version 2.0. des associations tontinières, avec le retour en force de la mutualisation des risques au sein de communautés soudées ?

Références

Joshua Angrist. The perils of peer effects. Labor Economics, 30, 98-108, 2014

Zakaria Babutsidze. Positions of French and German Banks in European interbank lending network. OFCE, Mars 2012.

Abhijit Banerjee, Arun Chandrasekhar, Esther Duflo & Matthew Jackson. Diffusion of Microfinance. Science, 341,

Smriti Bhagat, Moira Burke, Carlos Diuk, Ismail Onur Filiz & Sergey Edunov. Three and a half degrees of separation. Facebook Research, 2016.

Nicholas Christakis & James Fowler Social Network Sensors for Early Detection of Contagious Outbreaks. PLoS ONE. 5 (9): e12948. arXiv:1004.4792 2015

David Easley & Jon Kleinberg. Networks, Crowds, and Markets. Cambridge University Press. 2010.

Scott Feld. Why your friends have more friends than you do, American Journal of Sociology, 96 (6): 1464–1477, 1991.

Matthew Jackson. Social and Economic Networks. Princeton University Press, 2010.

Jure Leskovec & Eric Horvitz. Planetary-Scale Views on a Large Instant-Messaging Network. Microsoft Research, 2008.

Miller McPherson, Lynn Smith-Lovin & James Cook. Birds of a Feather: Homophily in Social Networks. Annual Review of Sociology. 27: 415–444, 2001.

James Moody. Race, School Integration, and Friendship Segregation in America. American Journal of Sociology, 107 (3): 679-716, 2001.

Wesley Perkins, Michael Haines & Richard Rice. Misperceiving the college drinking norm and related problems: a nationwide study of exposure to prevention information, perceived norms and student alcohol misuse. Journal of Studies on Alcohol 66 (4) : 470-478, 2005.

Ben Schiller. A Social Network For Insurance That Cuts Costs And Reduces Fraud. Fast Company, October 2013,

Brad Walker. How Peer-to-Peer Companies Are Transforming the Insurance Sector. The Street, Avril 2016,

Jason Wohlgemuth & Mihaela Matache. Small-World Properties of Facebook Group

Networks. Complex Systems, 23 (2014).

[i] Les données complètes sont en ligne sur https://snap.stanford.edu/data/email-Enron.html

[ii] https://www.friendsurance.com/ et https://www.inspeer.me/ respectivement