# On Wigner’s law (and the semi-circle)

There is something that I love about mathematics: sometimes, you discover – by chance – a law. It has always been there, it might have been well known by some people (specialized in some given field), but you did not know it. And then, you discover it, and you start wondering how comes you never heard about it before… I experienced that feeling this evening, while working on the syallbus for my course on copulas and extreme values. I discovered the so-called Wigner’s Semicircle Law (see e.g. Fan Zhang’s notes, or Fraydoun Rezakhanlou’s notes on that topic). Consider some $n\times n$ random matrice, with $n$ large (say 100) where elements are centered, such as a collection of random variable taking value $\{-1,1\}$ with equal probability. Then, eigenvalues can be visualized below

n=100
M=matrix(sample(c(-1,1),size=n*n,replace=TRUE),n,n)
E=eigen(M)$values plot(E,xlim=c(-11,11),ylim=c(-11,11)) Consider the symmetric matrix obtained from that matrix, $\frac{M+M'}{2}$ and more precisely, let us look at its eigenvalues, E=eigen(.5*(M+t(M)))$values

Then the distribution of those eigenvalues is the so-called semi-circle distribution

hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

Now, if we consider some $\mathcal{N}(0,1)$ distribution, instead of our binomial one, we got exactly the same

M=matrix(rnorm(n*n),n,n)
E=eigen(M)$values plot(E,xlim=c(-11,11),ylim=c(-11,11)) E=eigen(.5*(M+t(M)))$values
hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

Actually, it is a very general result, see the second chapter of an Introduction to Random Matrices by Greg Anderson, Alice Guionnet and Ofer Zeitouni, for instance. If entries of the random matrix are independent centred random variables, symmetric, such that higher moments exist, then this property is valid. That’s awesome, isn’t it? Because if the distribution has too heavy tails, then this property is no longer valid. For instance, if we consider a random matrix where entries have a Student distribution, we get something different…

M=matrix(rt(n*n,df=2.1),n,n)
M=M/sd(M)
E=eigen(M)$values plot(E,xlim=c(-11,11),ylim=c(-11,11)) E=eigen(.5*(M+t(M)))$values
hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

(here, I do normalize by the standard deviation to get something comparable with the previous graph, where variables were centered, with unit variance)

and if we consider a distribution with infinite variance, we get

M=matrix(rt(n*n,df=1.75),n,n)

we get

I guess I will get back on that property in my course!

# Random points on the Earth

The problem with puzzles is that you keep it in your head for days, until you find an answer. Or at least some ideas about a possible answer. This is what happened to me a few weeks ago, when a colleague of mine asked me the following question : Consider $n$ points uniformly distributed on a sphere. What is the probability that the $n$ points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?

Analogously, what is the probability to see the $n$ points on the Earth, at the same time, from somewhere in the galaxy ? (even extremely far away, so we can see a complete hemisphere) I wanted to use Monte Carlo simulations to estimate that probability, for some $n$. But it was difficult. I mean, given $n$ points on the sphere, in  can you easily determine if they lie on a common hemisphere, or not ? I did try with distance, or angle, but I could not find a simple answer. So I tried a technique I did learn a few years ago : if you cannot do something in dimension 3, try first in dimension 2.

Again, I could not find a simple answer. But there is a simple technique.

1. Draw $n$ points on the unit sphere, which simply means generate $n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$
2. Try to find $\theta\in[0,2\pi]$ such that, after a rotation (with angle $\theta$) all the points lie in the upper part (the North hemisphere, for instance)

So the question here is simply : is

$\max_{\theta\in[0,2\pi]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\geq 0)\right\}$

equal to $n$ ? Of course, from a computational point of view, it is slightly more complex, since this function is not differentiable.

n=5
Theta=runif(n)*2*pi
top=Vectorize(function(theta) sum(sin(Theta+theta)>=0))

So a simple strategy can be to compute those values on a finite grid, and to check if, for some $\theta$, all the points lie in the upper part

max(top(seq(0,2*pi,length=6001)))==n

Hence, with the following sample, all the points cannot be on a common hemisphere

set.seed(2)
Theta=runif(5)*2*pi

while, for another sample of points, it was possible

set.seed(7)
Theta=runif(5)*2*pi

With this simple code, we get get the probability, but only in dimension 2 (so far)

SIM=Vectorize(function(n) simul(n,1000))
plot(3:10,SIM(3:10))

In dimension 3, it is still possible to use also a polar representation. Things are easier to generate, but also it is simple to consider rotations. And again, a simple algorithm can be derived,

But there is a simple technique.

1. Draw $n$ points on the unit sphere, which simply means generate $2n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$ and $V_1,\ldots,V_n\sim\mathcal{U}([-\pi/2,\pi/2])$
2. Try to find $\theta$ and $\varphi$ such that, after a rotation (with angles $\theta$ and $\varphi$) all the points lie in some given part (say  $\{x\geq0\}$)

As mentioned by Dominique, using this technique, points are not uniformly distributed on the sphere. Instead we can use

1. Draw $n$ points on points from a trivaraite Gaussian distribution, and normalize it
2. Get the polar coordinates, and try to find $\theta$ and $\varphi$ such that, after a “rotation” (with angles $\theta$ and $\varphi$) all the points lie in some given part (say  $\{x\geq0\}$)

So the question here is simply : is

$\max_{(\theta,\varphi)\in[0,2\pi]\times[-\pi/2,\pi/2]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\cos(V_i+\varphi)\geq 0) \right\}$

(I am not sure about the set of angles for the rotations, so I tried a larger one, just in case). Again, it would be complex, or more complex than before, because we need here a joint grid. For instance

MZ=matrix(rnorm(n*3),n,3)
d=apply(MZ,1,function(z) sqrt(sum(z^2)))
X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d;
Theta=acos(Z)
Phi=acos(X/sqrt(X^2+Y^2))*(Y>=0)+(2*pi-acos(X/sqrt(X^2+Y^2)))*(Y<0)
top=function(theta,phi) sum(sin(Theta+theta)*cos(Phi+phi)>=0)
TOP=mapply(top,rep(seq(0,2*pi,length=1001),1001),rep(seq(-pi,pi,length=1001),each=1001))
max(TOP)==n

As we can see with the red curve, there might be some problems here. Because, (as mention in kmath327), this problem was solved in any dimension in Wendel (1962) with the following simple expression : in dimension $d$, the probability that the $n$ points (uniformly distributed on a sphere) lie on a same hemisphere is exactly

$p(d,n)= \frac{1}{2^{n-1}}\sum_{i=0}^{d-1} \binom{n-1}{i}$

p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))

Note that Leonard Savage proved (a few years before) that

$p(d,d+1)=1 -\frac{1}{2^d}$

(which can be obtained easily actually). I do not see what’s wrong with my Monte Carlo simulations… and if anyone has a nice Monte Carlo strategy to get that probability, I’d be glad to hear it !

# Probabilité et géométrie

Une des formules les plus importantes en probabilité (je trouve) est la “formule des probabilités totales” qui dit tout simplement que

$\mathbb{P}(Y=y)=\sum_x \mathbb{P}(Y=y,X=x)$

que l’ont peut aussi écrire, à l’aide de la formule de Bayes

$\mathbb{P}(Y=y)=\sum_x \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)$

Une des conséquences de ce résultat est la “law of total expectation“, souvent appelé théorème de double projection,

$\mathbb{E}(Y)=\sum_x \mathbb{E}(Y\vert X=x)\cdot\mathbb{P}(X=x)$

que l’on écrit souvent sous la forme raccourcie $\mathbb{E}(Y) =\mathbb{E}(\mathbb{E}(Y\vert X))$ (dans la formule de droite, le premier symbole est un espérance, c’est à dire une intégrale, et donc un nombre réel, alors que le second indique que l’on travaille sur une espérance conditionnelle, c’est à dire une variable aléatoire). Mais comme toujours dans les relations simplifiée, il faut savoir de quoi on parle.

La démonstration est simple

$\mathbb{E}(Y)=\sum_y y\cdot \mathbb{P}(Y=y)=\sum_y y \cdot \sum_x \mathbb{P}(Y=y, X=x)$

soit

$\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)=\sum_x\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)$

ou enfin

$\sum_x \mathbb{P}(X=x)\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)=\sum_x \mathbb{P}(X=x)\cdot \mathbb{E}(Y\vert X=x)$

tout simplement.

L’interprétation est encore plus simple (pour utiliser un exemple que j’aime de plus en plus) : si est le poids d’une personne $Y$, et $X$ le sexe, alors

$\mathbb{E}(Y)= \mathbb{E}(Y\vert X=\text{H})\cdot\mathbb{P}(X=\text{H})+\mathbb{E}(Y\vert X=\text{F})\cdot\mathbb{P}(X=\text{F})$

c’est à dire que le poids moyen d’une personne prise au hasard est un barycentre (damned, je fais déjà de la géométrie) entre le poids moyen des hommes et le poids moyen des femmes, les poids pour le calcul du barycentre étant liés aux proportions d’hommes et de femmes. Tout simplement.

Il existe une autre relation classique en probabilité (connue en statistiques sous le nom de “formule de décomposition de la variance”) qui dit que

$\text{Var}(Y)=\text{Var}(\mathbb{E}(Y\vert X))+\mathbb{E}(\text{Var}(Y\vert X))$

Là aussi, on a une écriture un peu simplificatrice, que l’on va essayer de décortiquer un peu. Pour ça, le plus simple est de faire un peu de géométrie. Car oui, dans les espaces de variables aléatoire, on peut faire de la géométrique. En particulier des projections (orthogonales). Mais avant, de faire des projections, il faut des distances, des angles, une notion d’orthogonalité, etc.

• Rappels de géométrie

Pour parler d’orthogonalité, il faut une notion de produit scalaire. Pour rappel, un produit scalaire, sur un espace $\mathcal E$, c’est défini par

• une notion de symétrie $<\vec x,\vec y>=<\vec y,\vec x>$
• de bilinéarité $<\alpha \vec x+\beta \vec y,\vec z>=\alpha<\vec x,\vec z>+\beta<\vec y,\vec z>$
• de positivité, $<\vec x,\vec x> \geq 0$
• et le produit scalaire est dit défini $<\vec x,\vec x>=0$ implique $\vec x=\vec 0$

De ce produit scalaire, on peut déduire une norme

$\|\vec x\|=\sqrt{<\vec x,\vec x>}$

En effet, on retrouve la propriété d’homogénéité, $\|\lambda\vec x\|=\vert\lambda\vert \|\vec x\|$ (c’est pour ça que l’on prend la racine carrée), et $\|\vec x\|=0$ si et seulement si $\vec x=\vec 0$. On a aussi l’inégalité de Cauchy-Schwarz,

$\|<\vec x,\vec y>\| \leq \|\vec x\|\|\vec y\|$

qui va impliquer l’inégalité triangulaire

$\|\vec x+y\|\leq \|\vec x\|+\|\vec y\|$

On dira que deux vecteurs sont orthogonaux, noté $\vec x\perp\vec y$ si $<\vec x,\vec y>=0$. Dans ce cas, on a la relation de Pythagore,

$\|\vec x+y\|^2= \|\vec x\|^2+\|\vec y\|^2$

Autre objet intéressant, la projection orthogonale. Rappelons tout d’abord que si $\mathcal F\subset \mathcal E$, on dira que $\vec x\perp\mathcal F$ si $\vec x\perp\vec y$y pour tout $\vec y\in\mathcal F$. Et on notera

$\mathcal F^\perp=\{\vec x\in\mathcal E;\vec x\perp\mathcal F\}$

On a alors des résultats intéressants. En particulier, pour tout $\vec y\in\mathcal E$, il existe un unique $\vec y_\star\in\mathcal F\subset \mathcal E$ tel que

$\|\vec y-\vec y_\star\|= \inf\{\|\vec y-\vec z\| ,\vec z\in\mathcal F\}$

On peut aussi montrer que $\vec y-\vec y_\star \in\mathcal F^\perp$. On parlera alors de projection orthogonale de $\vec y$ sur $\mathcal F$, et on pourra noter $\vec y_\star=\Pi_{\mathcal F}(\vec y)$. Tout cela est assez standard dans les espaces de dimension finie (on pensera à $\mathbb{R}^n$ pour avoir un peu d’intuition). On connait plein de choses sur les projections orthogonales, il suffit de penser dans $\mathbb{R}^n$ pour avoir un peu d’intuition. Par exemple, si $\mathcal{F}\subset\mathcal{G}$, alors

$\Pi_{\mathcal{F}}(\vec x)=\Pi_{\mathcal{F}}\big(\Pi_{\mathcal{G}}(\vec x)\big)$

Cette relation est vraie dans $\mathbb{R}^n$, mais aussi dans des espaces plus généraux. C’est ce qu’on appelle le théoreme de double projection (on qu’on devrait voir réapparaître sur les variables aléatoires).

• Géométrie dans les espaces de variables aléatoires

L’espace des variables aléatoires de variance finie – $L_2$ – peut être muni d’une telle opération, de produit scalaire,

$=\mathbb{E}(XY)$

La norme est alors $\|X\|=\sqrt{\mathbb{E}(X^2)}$. En fait, on voit qu’on s’avance sur un terrain glissant ici, car pour avoir un produit scalaire, il faudrait que $=0$ si et seulement si $X=0$, et pour avoir une norme il faudrait que $\|X\|=0$ si et seulement si $X=0$. Le soucis technique ici est que $\mathbb{E}(X^2)=0$ signifie que $\mathbb{P}(X=0)=1$ et pas $X(\omega)=0$ pour tout $\omega$. Bref, l’égalité est a comprendre au sens presque partout. Mais c’est un point de détail (ici). Techniquement, comme l’explique justement Williams (1991) dans le chapitre 6, $L_2$ est précisément l’espace $\mathcal L_2$ (classique en probabilité) quotienté par cette relation d’équivalence $\overset{p.s.}{=}$.

Un sous-espace de $L_2$ est l’espace des constante, $\mathbb{R}^n$, que je noterais aussi $s\{\boldsymbol{1}\}$. On notera que, comme rappelé dans tous les cours que je donne

$\mathbb{E}(Y) =\underset{c\in\mathbb{R}}{\text{argmin}}\{\mathbb{E}([Y-c]^2)\}$

c’est à dire que

$\mathbb{E}(Y)=\Pi_{s\{\boldsymbol{1}\}}(Y)$

L’espérance est la projection orthogonale sur l’espace des constantes. Et

$\text{Var}(Y) =\underset{c\in\mathbb{R}}{\text{min}}\{\mathbb{E}([Y-c]^2)\}$

Un autre sous-espace de $L_2$ qui sera intéressant est le suivant. Si $X\in L_2$, on notera

$s\{X\}=\{\tilde X=\psi(X)\in L_2,\psi:\mathbb{R}\rightarrow\mathbb{R}\}$

Par exemple, on peut considérer une variable de Bernoulli, $X=\boldsymbol{1}_A$$A\subset \mathbb{R}$. Dans ce cas, le sous-espace est équivalent au sous-espace des combinaisons linéaires,

$sl\{X\}=\{\tilde X=\beta_0+\beta_1 X;(\beta_0,\beta_1)\in \mathbb{R}^2\}\subset L_2$

Techniquement, je pense que

$\overline{\sigma\{X\}} \subset s\{X\}$

mais il semblerait qu’on ait égalité stricte… Mais comme je veux parler de projection sans parler de l’espace $\Omega$, on va continuer avec mon interprétation heuristique, et je ne parlerais pas de $\sigma\{X\}$. Je fais des projections sur des sous-espaces de $L_2$, pas sur des $\sigma$-algèbres.

Si $Y\in L_2$, on posera

$\mathbb{E}(Y\vert X)=\Pi_{s\{\boldsymbol{X}\}}(Y)$

qui sera la solution du problème de moindres carrés

$\mathbb{E}(Y\vert X) =\underset{X_\star\in s\{X\}}{\text{argmin}}\{\mathbb{E}([Y-X_\star]^2)\}$

On notera que $\mathbb{E}(Y\vert X) \in s\{X\}}$, et on retrouve ici l’écriture standard d’un modèle de régression, $\mathbb{E}(Y\vert X) =\psi(X)$. On parlera alors d’espérance conditionnelle. Qui est ici une variable aléatoire, par construction. En fait, l’unicité de la variable aléatoire est possible précisément parce que tout a l’heure, on définissait une égalité presque sure.

On avait rappelé le théoreme de double projection tout a l’heure. Et comme $s\{\boldsymbol{1}\}\subset s\{X\}$, alors

$\Pi_{s\{\boldsymbol{1}\}}(Y)=\Pi_{s\{\boldsymbol{1}\}}\big(\Pi_{s\{X\}}(Y)\big)$

ce qui se traduit par la relation bien connue

$\mathbb{E}(Y)= \mathbb{E}\big(\mathbb{E}(Y\vert X)\big)$

En économétrie (linéaire) et dans les premiers cours de séries temporelles, on introduit un autre type d’opérateur car on ne projette pas dans des espaces aussi gros que $s\{X\}$. On introduit souvent l’opérateur d’espérance linéaire, avec

$\text{EL}(Y\vert X)=\Pi_{sl\{{X}\}}(Y)$

où (comme on l’avait introduit tout à l’heure)

$sl\{X\}=\{\tilde X=\beta_0+\beta_1 X;(\beta_0,\beta_1)\in \mathbb{R}^2\}\subset L_2$

et pour les séries temporelles stationnaires,

$\text{EL}(X_t\vert X_{t-1},X_{t-2},\cdots)=\Pi_{\overline{sl\{\boldsymbol{X}\}}}(Y)$

Oui, si le processus est stationnaire, je peux régresser sur tout le passé, mais dans ce cas, je pense qu’il faut fermer l’espace afin de projeter dessus…

Bon, maintenant, on voit que mon interprétation géométrique est bancale… en effet, je ne peux définir mon espérance conditionnelle (voire mon espérance tout court) à l’aide de projections orthogonales que si je suis dans un espace muni qu’un tel opérateur. Et si $L_2$ est un espace de Hilbert, ce n’est pas le cas de $L_1$. Pourtant, nul besoin d’être dans $L_2$ pour définir une espérance conditionnelle… Je pense que l’astuce peut être de noter que si $X$ appartient à $L_1$ mais pas $L_2$, on peut quand même s’en sortir car $L_1$ est un espace de Banach et on a malgré tout des notions de convergences. Et $L_2$ est dense dans $L_1$. Bref, en faisant un peu de limite, pour peut étendre ce qu’on vient de faire pour les variables qui ne sont pas de carré intégrable…

Visuellement (il serait peut-être temps de faire des dessins, non ?) on a

On voit apparaître un paquet de relations sur ce dessin, que ce soit des doubles projections, mais aussi le théorème de Pythagore. Mais prenons notre temps…

• et la constante ?

Maintenant, si on regarde un peu détails, il y a des choses étranges. Par exemple $X\perp Y$ se traduit ici par $\mathbb{E}(XY)=0$. Mais ce n’est pas ce qu’on utilise classiquement comme notion d’orthogonalité. On a davantage l’habitude de voir

$\text{Cov}(X,Y)=\mathbb{E}(XY)-\mathbb{E}(X)\mathbb{E}(Y)=0$

Qu’est ce qu’on a raté ? L’idée est de poser $\vec x=X-\mathbb{E}(X)$, et $\vec y=Y-\mathbb{E}(Y)$, c’est à dire que l’on va travailler sur les variables centrées. On peut noter que

$\text{Var}(X)=\mathbb{E}([X-\mathbb{E}(X)]^2)=\| \vec x \|^2$

et on va définir la covariance entre X et Y a partir des vecteurs translatés,

$\text{Cov}(X,Y)=<\vec x,\vec y>$

et la corrélation comme

$\text{corr}(X,Y)=\frac{<\vec x,\vec y>}{\| \vec x \|\cdot \| \vec y \|} =\cos(\theta)$

où $\theta$ est l’angle entre les vecteurs $\vec x$ et $\vec y$. On dira que $X\perp Y$ si $\vec x \perp \vec y$ (au sens géométrique du terme), qui correspond a la relation classique $\text{corr}(X,Y)=0$). Le théroreme de Pythagore nous dit que si $\vec x \perp \vec y$, alors

$\| \vec x +\vec y \|^2=\| \vec x \|^2+\| \vec y \|^2$

$\text{Var}(X+Y)=\text{Var}(X)+\text{Var}(Y)$

si les variables sont orthogonales.

Maintenant, on était allé un peu plus loin tout a l’heure, avec des variables correspondant à des espérances conditionnelles. Par exemple, la formule de décomposition de la variance. On peut reprendre l’égalité du théorème de Pythagore, avec

$\vec x=\mathbb{E}(Y\vert X)-\mathbb{E}(Y)$

et

$\vec y=Y-\mathbb{E}(Y\vert X)$

On note facilement que $\vec x \perp \vec y$, car $\vec x\in s\{X\}$ alors que $\vec y\in s\{X\}^\perp$, par construction de la projection orthogonale. Aussi le théorème de Pythagore nous dit que

$\| Y-\mathbb{E}(Y) \|^2=\| \mathbb{E}(Y\vert X)-\mathbb{E}(Y) \|^2+\| Y-\mathbb{E}(Y\vert X) \|^2$

Pour le terme de gauche, c’est assez facile,

$\| Y-\mathbb{E}(Y) \|^2= \text{Var}(Y)$

Pour le premier terme de droite, la aussi, c’est facile, car $\mathbb{E}[\mathbb{E}(Y\vert X)]=\mathbb{E}(Y)$, par le théorème de double projection, et donc

$\| \mathbb{E}(Y\vert X)-\mathbb{E}(Y) \|^2=\text{Var}(\mathbb{E}(Y\vert X))$

Pour le dernier terme, a droite, c’est plus vicieux. Disons qu’on va l’identifier a une variance conditionnelle (je n’ai toujours pas défini formellement de variable aléatoire appelée “variance conditionnelle“), i.e.

$\| \vec Y-\mathbb{E}(Y\vert X) \|^2 = \mathbb{E}(\text{Var}(Y\vert X))$

On retrouve alors la formule de décomposition de la variance, en variance intra et variance inter.

Cette formule est classique en statistique, en régression (les résidus sont orthogonaux a $sl\{\boldsymbol{X}\}$, on aura donc une partie de la variance qui sera expliquée par nos variables explicatives $\boldsymbol{X}$ et une partie qui sera dite non expliquée) ou en crédibilité (je peux renvoyer cette fois à l’exemple des jeux de fléchettes de Philbrick (1982)). Promis, on reparlera de géométrie cet hiver, quand je ferais les lois multivariées en cours (avec cette fois des notions d’invariance, de rotations, de symétries, etc)

# Mathématiques de l’Assurance Non-Vie (2)

« Dans ce contexte d’incertitude, il est particulièrement réconfortant de revenir aux sources, aux fondamentaux, c’est-à-dire aux mathématiques et de rappeler que le risque naît de l’aléa et s’appréhende grâce aux développements les plus avancés du calcul des probabilités. » (Claude Bébéar, dans la préface du tome 1)

Après avoir été épuisé plusieurs semaines (voire plusieurs mois ?), une nouvelle impression du tome 2 de Mathématiques de l’Assurance Non-Vie, coécrit avec Michel Denuit, est de nouveau disponible. J’ai été surpris quand on nous a dit – en début d’année – que le second tome serait bientôt épuisé, et que plus d’un millier d’exemplaires avaient été vendus (ou donnés, ou volés… peu importe). Surpris qu’autant de monde pouvait être intéressé par un gros pavé (les deux tomes doivent faire plus de mille pages, au total) aussi sérieux. Comme le notait Gilles Bénéplanc dans la revue Risques, « s’ils sont riches et rigoureux (le titre présente l’avantage ne pas être trompeur quant au contenu), ces deux volumes d’actuariat se lisent toujours avec intérêt et même souvent avec plaisir ». Car effectivement, on a passé du temps à choisir le titre, et on est vite tombé d’accord sur le fait qu’il fallait insister sur le fait qu’il s’agissait de livre de Mathématiques (avant tout). Pour la petite anecdote, j’étais allé le chercher chez Gibert (juste après la publication du premier tome… après avoir fait un détour par les rayons de bandes dessinées, bien sûr). En vain. Là, un vendeur m’a expliqué que le livre était considéré comme traitant d’assurance, et que l’assurance était rangé dans les rayons de droit. Bref, notre livre vert était perdu au milieu des livres de Dalloz (alors que les livres de mathématiques financières se trouvaient dans la section mathématiques). Mais c’est un détail, tant qu’en le cherchant, on le trouve…

Je dois aussi dire que je suis très fier du livre. Très fier que des étudiants et des praticiens me disent le lire. Très fier d’apprendre que des étudiants au Brésil avaient photocopié (et traduit) l’introduction aux modèles de ruine (que l’on introduit avec l’approche de de Finetti, beaucoup plus intuitive je trouve…). Très fier quand des chercheurs que j’admire me disent avoir trouvé un petit résultat qu’ils cherchaient pour un papier de recherche (je pense a des discussions avec Alain, à Paris 1, sur la transformée d’Esscher). Je trouve aussi dangereux qu’autant de monde le lisent, car il est truffé d’erreurs. Mais après près de 6 mois de relecture complète, presque tous les jours, on fini par ne plus rien voir. Et malgré le travail incroyable de Christian, éditeur de la collection, il reste des coquilles (des coquillettes j’espère). Mais si on devait attendre que les livres soient parfaits avant de les publier, rien ne sortirait. Après presque dix ans, je pense que je réécrirais des paragraphes d’une autre manière, ou que je mettrais plus d’aspects algorithmiques. Mais ce qui est fait est fait, et j’en suis encore assez fier. OK, très fier.

Lorsque l’on a été contacté, on nous a demandé s’il fallait un  nouveau tirage, ou si une seconde édition (revue et corrigée) serait nécessaire. Soyons réalistes. Le livre n’est plus forcément à jour, et comme je l’ai dit, il reste beaucoup de coquilles, voire de fautes (il faudrait que l’on prenne le temps de les recenser). Mais une nouvelle édition ferait probablement doubler la taille du livre, et prendrait beaucoup de temps ! Sur les deux premiers chapitres, Michel a publié un livre complet sur les données de comptage en assurance. Sur les aspects computationnels, je vais éditer un livre à venir. Et s’il y a 10 ans, il y a avait encore peu de livres sur le sujet, plusieurs permettent maintenant de rentrer plus en détails sur certains aspects. Je pense au livre de Mario Wüthrich et Michael Merz sur le provisionnement (plus clair que le chapitre de notre livre, et plus complet). Pareil pour les chapitres de micro économie de l’assurance, ou sur les valeurs extrêmes ! On avait mis ensemble nos notes de cours, et on a pris plaisir à travailler ensemble sur ce projet. Et le livre reste pertinent pour nos cours. Il faudra juste aller voir ailleurs pour aller plus loin sur certains aspects ! Donc toutes nos excuses d’avoir juste demandé un retirage du livre, et de ne pas avoir pris le temps de faire une révision complète du livre…

# Newton-Raphson avec un dessin

Dans un précédant billet, publié il a quelques mois maintenant, je trouvais dommage que l’on oublie de citer des personnes qui ont introduit des concepts fondamentaux. En l’occurence, je pensais à Bruno de Finetti. En traînant sur Google Books l’autre jour, j’ai eu l’impression inverse, en trouvant dommage que l’on ne se souvienne que de la personne qui a introduit un concept, et pas celui qui en a compris la portée générale. Mais je devrais peut-être revenir un peu sur le contexte.

En travaillant avec un étudiant cet été, on devait résoudre des problèmes d’optimisation (en l’occurence l’estimation de paramètres dans un mélange de lois par maximum de vraisemblance) et je me souviens encore lui avoir dit, “c’est facile, tu programmes une descente de gradient, la méthode de Newton-Raphson, et c’est bon, on réglera l’histoire des contraintes plus tard” (oui, on avait un ensemble – simple – de contraintes qu’on pouvait intégrer en reparamétrant le programme d’optimisation). La méthode de Newton repose sur l’idée que pour trouver le zéro d’une fonction, on utilise une suite définie par récurence,

$a_{k+1} = a_k - \frac{f(a_k)}{f'(a_k)}$

La fameuse illustration de cette formule étant un dessin de la forme

De là, si on se souvient que la condition du premier ordre pour cherche un extremum d’une fonction revient à annuler la dérivée, on peut utiliser cette méthode pour chercher un extremum, avec une suite définie par

$a_{k+1} = a_k - \frac{f'(a_k)}{f''(a_k)}$

Cette méthode se généralise en dimension quelconque, avec la méthode dite de la descente de gradient

$\boldsymbol{a}_{k+1} = \boldsymbol{a}_k - [H f(\boldsymbol{a}_k)]^{-1} \nabla f(\boldsymbol{a}_k)$

Bref, ces méthodes sont très géométriques. Et j’avais toujours cru que les vieux mathématiciens étaient tous des géomètres. Aussi, j’espérais trouver plein de dessins dans les textes de Newton et Raphson. Mais non….

La méthode de Newton a été introduite par Isaac Newton dans De analysi per aequationes numero terminorum infinitas, écrit en 1669 et publié en 1711. Mais c’est surtout dans De metodis fluxionum et serierum infinitarum, écrit en 1671, traduit (du latin) et publié sous le titre Methods of Fluxions en 1736 que l’on va trouver l’idée de la méthode dite de Newton pour trouver les zéros d’une fonction à valeurs réelles. C’est ce texte que l’on peut trouver en ligne, numérisé. Mais Newton n’a aucunement proposé une méthode générale, car il ne lappliquait qu’aux seuls polynômes. A sa décharge, la notion de dérivée n’était pas (clairement) définie à cette époque… donc imaginer une descente en suivant la tangente, en 1710, ne serait peut-être pas raisonnable.

C’est Thomas Simpson qui proposa de généraliser cette méthode pour mettre en place une méthode de calcul itératif des solutions d’une équation non linéaire, en utilisant les dérivées (qu’il appelait fluxions, comme Newton).  En 1690,  Raphson publie une description de sa méthode dans Analysis aequationum universalis. Comme Newton, il va proposer une méthode de calcul récursif des approximations successives d’un zéro d’un polynôme (avec une petite nuance, car Newton essayait de construire une suite de polynômes, comme l’expliquent Ypma (1995) ou Deuflhard (2005))

Maintenant, si on regarde les publications, le dessin avec la tangente, on ne le voit pas ! Loin de là ! On voit des pages de calculs algébriques,

Il m’a fallu un certain temps pour découvrir un traité abordant ce problème avec des dessins. Pour les amateurs de mathématiques marseillaises, on peut découvrir un traité rafraîchissant, datant de 1768, publié d’ailleurs par un homme qui allait devenir par la suite maire de Marseille (à l’époque où les hommes politiques faisaient des maths, comme j’en parlais dans un vieux billet, avec une preuve du théorème de Pythagore proposé par un mathématicien amateur, qui allait devenir ensuite Président des États-Unis d’Amérique), monsieur Jean-Raymond Pierre  Mouraille.

Ce traité est passionnant, avec un style d’une autre époque

Mourailles se positionne clairement dans la lignée de Newton, avec quelques bémols toutefois,

Le plus intéressant, ce sont les pages d’annexes, remplies de dessins, où enfin, on peut découvrir le dessin avec la suite définie par récurrence, et la descente le long de la tangente,

Bref, tout le monde pense à ces dessins quand il entend parler de la méthode de Newton-Raphson, même si aucune représentation géométrique n’est présentée.

Je ne sais pas si Jean-Raymond Pierre  Mouraille a le premier a avoir représenté la recherche de zéros d’une fonction de cette manière, mais il est dommage qu’on n’en parle pas plus souvent…

# Generating a Markov chain vs. computing the transition matrix

A couple of days ago, we had a quick chat on Karl Broman‘s blog, about snakes and ladders (see http://kbroman.wordpress.com/…) with Karl and Corey (see http://bayesianbiologist.com/….), and the use of Markov Chain. I do believe that this application is truly awesome: the example is understandable by anyone, and computations (almost any kind, from what we’ve tried) are easy to perform. At the same time, some French students asked me specific details regarding some old lectures notes on Markov chains, and on some introductory example I used as a possible motivation: the stepping stone algorithm. In the notes, I just mentioned the idea of this popular generic algorithm (introduced in Sawyer (1976)) and I use simulations to show – visually – how it works. Again, it was just to motivate the course which actually did focus on the theory of Markov Chains. But those student wanted more, like how did I get the transition matrix, for instance. And that is actually not a simple question, from a computational perspective. I mean, I can easily generate this Markov Chain, but writing explicitly the transition, that was another story. Which took me a bit longer. In a very specific case…

But let us get back to the roots, and to the stepping stone algorithm. At least, one of them (the one I used in my notes) because it looks like there are several algorithm. We do consider a grid, say $h\times h$, with some colors inside, say $k$ possible colors. Each cell of the grid has a given color. Then, at some stage, we select randomly one cell in the grid, and it will take the color of one of its neighbor (some kind of absorption, or mutation). This is, more or less, what is also detailed in some lecture notes by James Propp (see also e Sato (1983) or Zähle et al. (2005) for more theoretical details about that Markov chain). This is extremely simple to generate (that’s what I did in my notes, with very big grids, and a lot of colors). But what if we want to write the transition matrix ?

First of all, we need to define the state space. Basically, we do have $h^2$ cells, each of them has one color, chosen among $k$. Which gives us $k^{h^2}$ possible states…. And that can be large. I mean, if we consider the smallest possible grid (that might be interesting), say $h=3$, and only $2$ colors, then we talk about $2^9=512$possible states. That is large, not huge. But we should keep in mind that we have to compute a transition matrix, that would be a matrix with $512\times 512=262,144$ elements. More generally, we talk about writing down matrices with $k^{2\cdot h^2}$ elements. If we want black and white $4\times 4$ grids, that would mean a matrix with $2^{2\cdot 4^2}=2^{32}=4,294,967,296$ which mean 4 billion elements ! And if we consider an red-green-blue $3\times 3$ grid, we have to explicit a matrix with $387,420,489$ i.e almost 400 million elements. So, let’s face it: we can only work with $3\times 3$ bi-color grids.

So let’s try… The good thing is that it can be related to work I’ve been doing recently on binomial recombining trees (binomial being related to bi-color). First of all, our grid will be describes as follows

> h=3
> M=matrix(1:(h^2),h,h)
> M
[,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

with two colors

> color=c("red","blue")

Then, we should look for neighbors, or derive an neighborhood matrix,

> d=function(i,j) dist(rbind(c((i-1)%/%h,(i-1)%%h),
+                            c((j-1)%/%h,(j-1)%%h)))
> Neighb=matrix(Vectorize(d)(rep(1:(h^2),each=h^2),
+                            rep(1:(h^2),h^2)),h^2,h^2)
> trunc(Neighb*100)/100
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 1.00 2.00 1.00 1.41 2.23 2.00 2.23 2.82
[2,] 1.00 0.00 1.00 1.41 1.00 1.41 2.23 2.00 2.23
[3,] 2.00 1.00 0.00 2.23 1.41 1.00 2.82 2.23 2.00
[4,] 1.00 1.41 2.23 0.00 1.00 2.00 1.00 1.41 2.23
[5,] 1.41 1.00 1.41 1.00 0.00 1.00 1.41 1.00 1.41
[6,] 2.23 1.41 1.00 2.00 1.00 0.00 2.23 1.41 1.00
[7,] 2.00 2.23 2.82 1.00 1.41 2.23 0.00 1.00 2.00
[8,] 2.23 2.00 2.23 1.41 1.00 1.41 1.00 0.00 1.00
[9,] 2.82 2.23 2.00 2.23 1.41 1.00 2.00 1.00 0.00
> Neighb=(Neighb<2)&(Neighb>0)
> Neighb
[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
[1,] FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE
[2,]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE
[3,] FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE
[4,]  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
[5,]  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
[6,] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE
[7,] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE
[8,] FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
[9,] FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE

Now, let us explicit our 512 possible states.

> n=h^2
> states=function(x){
+   Base.b=rep(0,n)
+   ndigits=(floor(logb(x,base=length(color)))+1)
+   for(i in 1:ndigits){
+     Base.b[n-i+1]=(x%%length(color))
+     x=(x %/% length(color))}
+   return(Base.b)}
> M=Vectorize(states)(1:(length(color)^n-1))
> liststates=data.frame(rbind(rep(0,h^2),t(M)))
X1 X2 X3 X4 X5 X6 X7 X8 X9
1  0  0  0  0  0  0  0  0  0
2  0  0  0  0  0  0  0  0  1
3  0  0  0  0  0  0  0  1  0
4  0  0  0  0  0  0  0  1  1
5  0  0  0  0  0  0  1  0  0
6  0  0  0  0  0  0  1  0  1

(for the first six, with 0/1 digits instead of colors). For instance, if we look at a specific one, it is possible to plot the grid, using

> plotsteps=function(u){
+   plot(0:h,0:h,col="white",xlab="",ylab="",axes=FALSE)
+   for(i in 0:(h^2-1)){
+   x=i%/%h
+   y=i%%h
+   polygon(x+c(1,.1,.1,1),y+c(1,1,.1,.1),
+   col=color[as.numeric(u)[i+1] + 1])
+   text(x+.45,y+.45,i)
+   }}

Here,

> plotsteps(liststates[100,])

Then, given one state, let us see what could happen next,

• let us compute all connected states: all states where we can end up in if we change one cell
• we have to check, for each connect state which cell did change
• we should compute probabilities to reach those 9 states, based on the fact that each of the cell is chosen with the same probability, and the fact that probability to change the color is based on the colors around.
• if some states cannot be reached (if a cell is surrounded by elements of the same color, so it cannot change its color), then, we should remove then from the list of reachable (possible) states.

The code will be something like the following

> listneighbour=function(i){
+   start=liststates[i,]
+   difference2only=function(j) {
+     w=which(liststates[j,]!=liststates[i,])
+     return((length(w)==1))}
+   possible=which( Vectorize(difference2only)(1:nrow(liststates))==TRUE )
+   P=function(j){
+     L=liststates[i,which(Neighb[which(liststates[j,]!=liststates[i,]),]==TRUE)]
+     T=table(as.numeric(L))
+     T=T[as.character(0:(length(color)-1))]
+     T[is.na(T)]=0
+     return(as.numeric(T)/sum(T))
+   }
+   probability=Vectorize(P)(possible)
+   W=NULL
+   for(j in possible) W=c(W,which(liststates[j,]!=liststates[i,]))
+   I=1-liststates[i,W]+1
+   vp=diag(probability[as.numeric(I),])
+   vproba=0*vp
+   if(sum(vp)!=0) vproba=vp/sum(vp)
+   return(list(
+     color=liststates[i,W],
+     absorb=W,
+     possible=possible,
+     probability=probability,
+     prob=vproba))
+ }

For instance, if we start from state 100 (here, on the right)

> listneighbour(100)
$color X3 X4 X8 X9 X7 X6 X5 X2 X1 100 1 1 1 1 0 0 0 0 0$absorb
[1] 3 4 8 9 7 6 5 2 1

$possible [1] 36 68 98 99 104 108 116 228 356$probability
[,1] [,2] [,3]   [,4]   [,5] [,6] [,7] [,8]   [,9]
[1,]    1  0.8  0.6 0.6667 0.3333  0.4  0.5  0.6 0.6667
[2,]    0  0.2  0.4 0.3333 0.6667  0.6  0.5  0.4 0.3333

$prob [1] 0.17964072 0.14371257 0.10778443 0.11976048 0.11976048 [6] 0.10778443 0.08982036 0.07185629 0.05988024 Let us look more specificaly at the 99th state (which appears above as a state that could be reached from the 100th), > liststates[99,] X1 X2 X3 X4 X5 X6 X7 X8 X9 99 0 0 1 1 0 0 0 1 0 If we plot it (here on the right, again), we get > plotsteps(liststates[99,]) Clearly, here, the cell in the upper corner (number 9) changed from blue to red. Now, about the probability… The probability to select cell 9 is 1/9, and given that cell 9 is chosen, the probability to go from blue to red is 2/3 (the cell is surrounded by 2 red cells, and 1 blue cell). The probability to remain blue is then 1/3. Those are the probabilities computed by our function (the table with two rows, one per color). In order to get a better understanding on the meaning of the last line, with some sort of probabilities), let us look at the following (simpler) example. > liststates[2,] X1 X2 X3 X4 X5 X6 X7 X8 X9 2 0 0 0 0 0 0 0 0 1 that can be visualized on the right (on the right). Here, > listneighbour(2)$color
X9 X8 X7 X6 X5 X4 X3 X2 X1
2  1  0  0  0  0  0  0  0  0

$absorb [1] 9 8 7 6 5 4 3 2 1$possible
[1]   1   4   6  10  18  34  66 130 258

$probability [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1 0.8 1 0.8 0.875 1 1 1 1 [2,] 0 0.2 0 0.2 0.125 0 0 0 0$prob
[1] 0.65573770 0.13114754 0.00000000 0.13114754 0.08196721
[6] 0.00000000 0.00000000 0.00000000 0.00000000

Things are pretty simple here

• if we chose cells $\{1,2,3,4,7\}$, then nothing change, since all the neighbors have the same color. So if we want to focus on changes (or say run the algorithm until the first color change, then choosing those cells is a waste of time)
• if we chose cells $\{5,6,8\}$, then it could be possible to change the color. And actually, $\{5\}$ is different from $\{6,8\}$ (since it does have much more neighbors)
• if we chose cell $\{9\}$, then definitively, the color will change, since all neighbors have the other color here,

Now, the probability to select cell $k$ given that there was a color change would be, if $k$ is in $\{9\}$

$\mathbb{P}(k)\propto \frac{3}{3}=1$

while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so

$\mathbb{P}(k)\propto \frac{1}{5}$and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so

$\mathbb{P}(k)\propto \frac{1}{8}$

And for the other, $\mathbb{P}(k)\propto 0$. So, it comes – since we assume that cells are drawn independently, and with the same probability, if $k$ is in $\{9\}$

$\mathbb{P}(k)= \frac{1 \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{40}{61}$

while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so

$\mathbb{P}(k)= \frac{\frac{1}{5} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{8}{61}$

and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so

$\mathbb{P}(k)= \frac{\frac{1}{8} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{5}{61}$

Which are exactly the probability computed above. The point is that we compute probabilities given that a color change will actually occur. The good point is that it should faster convergence to some limiting distribution. If any.

What about our transition matrix ? Well, using a simply loop, we should get it easily

> M=matrix(0,nrow(liststates),nrow(liststates))
+ for(i in 1:nrow(liststates)){
+ L=listneighbour(i)
+ if(sum(L$prob)!=0){ + j=L$possible
+ M[i,j]=L$prob + } + if(sum(L$prob)==0){
+ j=i
+ M[i,j]=1
+ }
+ }

One can check that this matrix satisfies some properties of transition matrices. For instance, the sum per row is one,

> sum(apply(M,1,sum)!=1)
[1]  0

Remember that this matrix is big, so I will not print if here. But trust me, it works (it might take a while on an old laptop, but anyone can do it). Now, if we want to visualize some paths of that chain, we can use the following algorithm. First, we need a starting point, that can be chosen randomly,

> j=sample(1:nrow(liststates),size=1)

or using a given colored grid, say

> j=100

Then we plot it,

> plotsteps(liststates[j,])

Now, the code within the loop is here

> d=rep(0,nrow(liststates))
> d[j]=1
> d=d%*%M
> j=sample(1:nrow(M),size=1,prob=d)
> plotsteps(liststates[j,])

Here are some examples. And indeed, we end up either with all cells in blue, or all cells in red.

Now, do we have to compute that transition matrix to produce those graph (and to generate that Markov chain) ? No. Of course not… At each step, I use a Dirac measure, and use the transition matrix just to get the probability to generate then the next state. Actually, one can write a faster and more intuitive code to generate the same chain… But I should probably keep that for another post…

I will give a talk tomorrow morning at the Mathematical  Finance Days, organized in HEC Montréal Monday and Tuesday, on Advanced methods in trees with (as mentioned in the subtitle of the first slide) a some thoughts on teaching mathematical finance. It is mainly a survey on advanced tools, based on the idea expressed in Price (1996),

The paper that showed that European option pricing could be put on a rational mathematical basis was Black and Scholes published in 1973. It was so revolutionary that the authors had to submit it to a number of journals before it was accepted. Although there are now numerous approaches to the result, they mostly require specialized methods, including Ito calculus and partial differential  equations, and perhaps Girsanov theory and Feynman-Kac methods. But it is the binomial method due initially to Sharpe and substantially extended by Cox, Ross, and Rubinstein that made the theory of option pricing accessible to everyone with limited mathematical background.  Even though it requires only routine algebraic manipulations, the method is still able to elucidate many of the ideas behind the full theory. Furthermore, all the surprising results mentioned in the opening can be located in this approach. For these reasons it is usually the first method presented in text books and finance courses; we shall follow this trend and step through it. The binomial method is, however, much more than a pedagogical breakthrough, since it allows for the development of numerical approximation methods for a wide range of options for which there are no  known analytic solutions.

Some recent results, obtained in work in progress with colleagues in combinatorial analysis will also be mentioned at the end of the talk (slides can be downloaded in a pdf format, with animations)

I will also be chairing the Numerical Methods session.

# Pavage et magie, rotations et symétries

(oui, mon fils voulait une rime dans le titre… promis on va tenter de faire une comédie musicale à partir de mes billets). La semaine passée, un article intéressant était publié (http://www.bbc.co.uk/news/…) sur sciences et magie. Cet article m’a fait penser à une réflexion que je m’étais faite lors de notre visite au MoMaths, cet hiver. Pour revenir sur l’histoire que j’avais raconté il y a quelques mois, sur la fin de la visite, les enfants avaient joué avec les pavages, à la Escher (les pièces étaient magnétiques, et s’assemblaient sur un grand tableau blanc)

Le soucis c’est qu’à aucun moment je n’ai pu voir d’explication du  pourquoi ça marche ? ou pour être plus explicite comment a été construit le singe pour que tout s’emboîte aussi bien. Car il y a un truc… La magie ça n’existe pas ! ce sont juste des maths (c’est un peu ce que racontent Persi Diaconis et Ron Graham dans Magical Mathematics, même si j’avoue ne pas l’avoir fini, faute de temps, et d’avoir laissé des chapitres en attente). Mon fils travaillant les rotations et les symétries, à l’école, j’ai voulu qu’on prenne le temps de revenir sur les pavages. Mais avant de me lancer avec lui, j’ai essayé de m’entraîner un peu tout seul (dimanche après midi, pendant le temps calme).

La difficulté (ou la magie) est qu’il y a des symétries dans tous les sens. Regardons un exemple simple (ou disons pas trop compliqué) avec des chinois (identifiés par les chapeaux dits chinois). En regardant rapidement, on note que le triangle ci-dessous est reproduit à l’infini, par un jeu de rotations. Il suffit donc d’avoir un tel triangle. Le soucis est que le découpage (parfaitement symétrique) n’est pas trivial (il est délicat de faire un personnage en faisant un découpage dans un triangle, même si c’est la méthode retenue dans l’activité mentionnée sur https://irem.univ-lille1.fr/…)

Le secret part de l’observation suivante: sur le graphique de gauche, en dessous, les trois points (qui sont des intersections des trois personnages, ou disons des trois couleurs) sont les sommets d’un triangle équilatéral. C’est là que la magie va opérer, car c’est cette figure (avec une autre, juxtaposée) que l’on va tenter de reproduire. En fait, un personnage (disons un petit chinois bleu) c’est deux triangles collés (on doit pouvoir parler de parallélogramme isocèle).

On ne le voit pas très bien, mais pourtant, au départ, on avait un parallélogramme isocèle (dont les contours sont tracés en noir) qui était  tout bleu. Bien entendu, on a découpé des parties et en les a envoyé ailleurs (par rotation autour d’un sommet). Le plus simple est de construire la figure, à la main. On part donc du parallélogramme suivant, tout bleu :

Oui, j’ai repris les fonctions de base utilisées l’autre jour pour expliquer R aux enfants,

source("http://freakonometrics.free.fr/RforKIDZ.R")

En fait, j’ai souhaité définir une fonction créant un polynôme à partir d’une liste de points (ou de sommets).

polyg=function(L,couleur="blanc",tour=couleur,ep=1,st=1,...){
b=unlist(L)
n=length(b)
polygon(b[seq(1,n,by=2)],b[seq(2,n,by=2)],col=couleurc(couleur),border=couleurc(tour),lwd=ep,lty=st,...)}

La fonction pour faire une rotation d’une liste de points est ici

rotation=function(L,angle=90,centre=L[[1]]){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]-centre[1],b[seq(2,n,by=2)]-centre[2])
Y=matrix(c(cos(angle/180*pi),sin(angle/180*pi),-sin(angle/180*pi),cos(angle/180*pi)),2,2)%*%t(X)
Z=Y+centre
L2=L
for(i in 1:(n/2)){L2[[i]]=Z[,i]}
return(L2)}

On a alors ici les quatre sommets

A=c(4,0)
B=c(4,10)
C=c(4+sqrt(3)*5,5)
D=c(4+sqrt(3)*5,15)

et on peut tracer le parallélogramme avec

L1=list(A,B,D,C)
polyg(L1,"bleu",tour="noir",ep=3)

On peut alors jouer à découper des parties, puis à les pivoter (on va faire des rotations, quoi). Par exemple, le prend le triangle en haut à droite, et je pivote de 240 degrés autour d’un des sommets du parallélogramme. C’est la figure ci-dessus à droite.  Le triangle, il est défini par le code suivant

E=c(10,10)
F=c(7,10+sqrt(3))
L1=list(B,E,F)
polyg(L1,tour="rouge",ep=3)

ce qui enlève la portion – en traçant un triangle blanc – et pour faire la rotation, j’utilise le code suivant

L2=rotation(L1,240)
polyg(L2,"bleu","rouge",ep=3)

(la seconde ligne permet de tracer le triangle bleu). Il va de soi qu’on peut découper des objets plus complexes que des triangles, voire non connexes. Par exemple, je peux découper la forme verte ci-dessous, puis la pivoter par rapport au sommet en haut à gauche, de 240 degrés.

G=c(10,2*sqrt(3)+10)
H=c(4+sqrt(3)*5,2*sqrt(3)+10)
F1=c(8,10+2/sqrt(3))
F2=c(8,10+4/sqrt(3))
L1=list(B,E,F1,F2,B,G,H,D,G)
polyg(L1,"vert",tour=NA)
L2=rotation(L1,240)
polygL(L2,"bleu",tour=NA)

(il n’est pas nécessaire de revenir au point B, mais peu importe). Ou la figure jaune, qui va pivoter autour du sommet en bas à droite, de 120 degrés.

F3=c(8,10-2/sqrt(3))
F4=c(8,10-4/sqrt(3))
G2=c(10,-2*sqrt(3)+10)
H2=c(4+sqrt(3)*5,-2*sqrt(3)+10)
L1=list(H,G,F2,F1,E,F3,F4,G2,H2,H)
polyg(L1,"jaune",tour=NA)
L2=rotation(L1,120,centre=C)
polygL(L2,"bleu",tour=NA)

Avec des ciseaux, ça sera plus simple, mais le but est aussi de voir que programmer, c’est facile. En gros, on a deux fonctions de base: une qui trace un polynôme et l’autre qui fait des rotations.

Si on fait les deux en même temps (à peu près), on enlève toute la partie en haut à droite, et on la répartie: une partie à gauche, et une partie en bas à droite.  On a alors une forme qui ressemble à un chinois. Non ? On peut aller creuser un peu sous les bras, pour améliorer les mains, mais on y est presque, me semble-t-il…

On peut continuer à améliorer, par exemple en enlevant des parties, puis en faisant des rotations par rapport à un des sommets initiaux. Comme le petit triangle jaune, qui va pivoter de 120 degrés par rapport au sommet en bas à droite. Puis le triangle vert par rapport au sommet en bas à gauche (on note qu’on pivote par rapport aux sommets du triangle équilatéral qu’on avait repéré, initialement)

S=c(8-.4,-2*sqrt(3)+10-.4)
L1=list(F4,G2,S,F4)
polyg(L1,"jaune")
L2=rotation(L1,120,centre=C)
polyg(L2,"bleu")

Bon, là on fait une pause, parce qu’on a, par morceaux, notre chinois ! On un paquet de lettre (voire un paquet de sommets qu’il va falloir nommer, suite à nos rotations). On les nomme avec un code qui ressemble à (oui, on aurait pu pu les nommer dès le départ, lors des rotations)

rE=rotation(list(E),240,centre=B)[[1]]
rF1=rotation(list(F1),240,centre=B)[[1]]
rH=rotation(list(H),240,centre=B)[[1]]

Une fois qu’on a identifié nos sommets, on a notre chinois,

chinois=list(B,rE,rF1,LG,S2,LF2,rH,A,r2H,r2G,r2S2,
r2F2,r2F1,r2E,rF3,rF4,r2S,rG2,rH2,C,H2,G2,S,F4,F3,E,B)

Bon, ensuite, reste à paver, ce qui est est le plus simple. Par exemple, on peut faire des rotations,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")

voire définir une translation,

translation=function(L,h){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]+h[1],b[seq(2,n,by=2)]++h[2])
L2=L
for(i in 1:(n/2)){L2[[i]]=X[i,]}
return(L2)}

puis translater notre chinois,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")
polyg(translation(r1chinois,15/sqrt(3),15),couleur="jaune")

Il manque l’autre, le vert. On peut vérifier qu’il se superpose parfaitement :

Bon, ben ça y est, on a constitué notre pavage. Reste à travailler un peu le code, essentiellement sur la forme, et on pourra jouer avec les enfants, à faire les lézards, des singes, et toutes sortes de pavages possibles, à la Escher. En notant qu’on peut utiliser d’autres formes de bases plus complexes, comme expliqué sur http://xavier.hubaut.info/…

# Les experts

Il y a un an, presque jour pour jour, je finissais un rapport que je remettais pour un litige qui devait passer devant un tribunal. J’avais été sollicité pour faire des simulations de scénarios, et étudier des lois de temps d’arrêt et des problèmes de probabilités de ruine. Bref, j’étais un peu devenu un expert. Peut-être pas exactement un expert judiciaire au sens strict, car j’étais sollicité par une des parties, et non pas pas un juge (est-ce toujours le cas ?). Et je n’ai pas non plus prêté (solennellement) le serment “Je jure d’apporter mon concours à la Justice, d’accomplir ma mission, de faire mon rapport, et de donner mon avis en mon honneur et en ma conscience”. Même si c’est ainsi que j’ai mené ma mission. Mais peu importe (pour des raisons de confidentialité, je n’en dirais pas plus sur cette mission). C’est toujours valorisant de se dire que les choses théoriques que l’on enseigne peuvent servir dans la vraie vie !

Hier, au détour d’une conversation avec mart1oeil, goulu ‏et tomroud, j’ai découvert zythom et son blog, http://zythom.blogspot.ca/, qui est le “blog d’un informaticien expert judiciaire“, comme l’indique l’entête. Et hier, zythom posait la question un expert judiciaire peut-il être mauvais dans votre domaine ? Je reviens sur son billet car le billet commence par une recherche de définition de ce que peut être un expert. Et comme toujours la page de Wikipédia consacrée à l’expert propose une lumière intéressante. “L’expert n’est pas simplement celui qui sait, sur un champ délimité de savoir. Son expérience reconnue lui permet d’apporter une réponse argumentée à une demande d’expertise. Il faut le différencier du savant et aussi du spécialiste.” Et zythom prend alors l’image du médecin de famille ou du médecine généraliste, en opposition au spécialiste. Bref, tout d’un coup, hier un monde s’est effondré ! Je suis probablement pas un expert !

Cela dit, j’aurais pu m’en douter plus tôt. En fait, je crois que je le savais depuis des années, mais que je me voilais la face.

Il y a quelques années maintenant, j’enseignais un cours de méthodes numériques en finance (j’avais même tapé des notes de cours). A un moment, on montre que le prix $V$ vérifie une belle équation, qui permettra d’arriver à la formule de Black-Scholes, qui permet d’avoir le prix d’une option dite européenne. L’équation est de la forme

$\displaystyle \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0$

En bricolant un peu (avec quelques changements de variables), on peut se ramener une équation de la forme

$\frac{\partial u(t,x)}{\partial t}-\kappa \frac{\partial ^2 u(t,x)}{\partial x^2}=f(t,x)$

Tous ceux qui ont fait un peu de maths reconnaîtrons une équation classique (datant de plus de 200 ans) appelée équation de la chaleur (les plus curieux pourront aller voir les notes de cours de Nicolas Vandewalle sur le sujet, entre autres). Et dans mon cours, je montrais comment résoudre numériquement cette équation, avec différents schémas, comme un schéma explicite

$\frac{u_{i+1,j}-u_{i,j}}{\Delta t}-\kappa \frac{u_{i,j+1}-2 u_{i,j}+u_{i,j-1}}{h^2}=f_{i,j},$

en insistant sur les conditions numériques pour que la résolution n’explose pas (comme sur le dessin ci-dessous)

Il fallait pour cela que le ratio

$\lambda=\kappa\frac{\Delta t}{h^2}$

soit inférieur à 1/2. Etc.

Bref, je suis capable de résoudre numériquement l’équation de la chaleur depuis des années, voire d’enseigner comment le faire. N’empêche que je lorsque je réchauffe mon lunch au micro-onde, tous les midis, je me brule, en ayant en même temps un morceau qui est encore congelé dans le plat ! C’est pénible ! Vivement que je devienne un expert !

# Martingale et journalisme scientifique

Être journaliste scientifique ne doit pas être facile. J’imagine qu’il faut être à l’écoute des nouvelles scientifiques, et d’informer, aussi justement que possible. C’est ce que fait avec brio Pierre Barthélémy (aka @PasseurSciences sur Twitter) dans sa chronique hebdomadaire sélection scientifique de la semaine. Il essaye ainsi de parler de sciences dans un journal qui a (trop) souvent confondu science et technologie (technologie étant aussi souvent un mot savant utilisé pour masquer de la publicité pour des appareils technologiques avancés). Et comme l’espace qui lui est imparti est restreint (c’est le moins qu’on puisse dire), la plupart des informations passent par son blog http://passeurdesciences.blog.lemonde.fr/, qui est riche ! Sylvestre Huet fait aussi ce travail pour Libération, avec son blog http://sciences.blogs.liberation.fr/ (là encore, on peut lire bien plus de choses en ligne que que dans la version papier, ce qui est d’autant plus intéressant pour les lecteurs de l’autre coté de l’océan atlantique). Sylvestre Huet arrive à parler géologie, puis démographie, avant d’évoquer dans un article truffé de références les réformes des instances qui régentent la recherche en France. Oui, les deux blogs sont malgré tout très français. Pour une vision plus internationale, on pourra citer l’admirable travail de Pascal Lapointe, par exemple (aka @paslap) sur http://sciencepresse.qc.ca/.

Mais si je commence à mettre les pieds dans le plat, je poserais la question de la légitimité d’un journaliste scientifique, à l’heure où la blogosphère scientifique explose. Par exemple, si je veux apprendre des choses en physique, ou en biologie, je sais que le blog de Tom Roud (que l’on peut suivre sur Twitter @tomroud) http://tomroud.cafe-sciences.org est une source infinie d’information (ok, peut-être pas “infinie” car comme nous tous, il a une vie en dehors de son blog, mais disons qu’il y a de quoi lire). En sciences humaines, il existe des centaines de blogs hébergés sur http://hypotheses.org/ tenus par des universitaires, ou http://www.cafe-sciences.org/ pour des blogs de sciences, en français. Les bloggers scientifiques s’expriment – le plus souvent – en restant dans leur champ d’expertise. Et c’est tant mieux. Le journaliste scientifique, lui, se doit de parler de tous les sujets scientifiques. Et c’est là que l’exercice devient délicat, car le journaliste se doit d’être critique. Compte tenu de la concurrence féroce qui existe dans le monde académique, on nous pousse à faire croire qu’on vient d’écrire l’article qui va révolutionner la science, sinon le monde. Que notre approche est novatrice, et qu’en plus, on vient de montrer ce qui pourrait ce qui pourrait être le sain Graal dans notre communauté. C’est le travail du journaliste de se demander si c’est vrai….

J’en arrive au point de mon billet. Le dernier article  le hasard, martingale boursière ? (en ligne sur http://lemonde.fr/sciences/…) de Pierre Barthélémy m’a un peu agacé. Pas sur le thème abordé, mais l’impression générale qui m’est restée après avoir lu l’article (samedi matin). Pour l’histoire complète, c’est ma femme qui est tombée la première sur l’article, et qui a été étonnée de lire un article pareil dans Le Monde, en tant que mathématicienne (et probabiliste). Et quand je l’ai lu, c’est en tant qu’économiste que je suis resté sans voix. En résumé (trop succinct, j’en convient), on nous résume l’article de physiciens italiens qui “expliquent qu’au grand dam des traders et autres analystes financiers, en quête perpétuelle de justifications rationnelles aux fluctuations des cours boursiers, les marchés demeurent obstinément imprévisibles“. Pour parcourir l’article mentionné, mais non cité – probablement car il s’agit de la version papier ? – l’article est en ligne sur http://arxiv.org/1303.4351  (trouvé via http://improbable.com/…). On retrouve cette idée dans le texte sous la forme suivante “our main result, which is independent of the market considered, is that standard trading strategies and their algorithms, based on the past history of the time series, although have occasionally the chance to be successful inside small temporal windows, on a large temporal scale, perform on average not better than the purely random strategy, which, on the other hand, is also much less volatile. In this respect, for the individual trader, a purely random strategy represents a costless alternative to expensive professional financial consulting, being at the same time also much less risky, if compared to the other trading strategies.” On apprend dans l’article paru dans Le Monde que “pour Pluchino et compagnie, le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre“. Damned, rien que ça ? Pour les amateurs de lyrisme (il doit y en avoir dans les lecteurs du monde, moins parmi les lecteurs plus scientifiques qui aiment les phrases moins alambiquées), je passe les moments où Pierre Barthélémy est clairement plus journaliste que scientifique, comme lorsqu’il nous explique que ces quatre chercheurs “viennent d ‘attaquer à l’acide de l’aléatoire” les modèles financiers. Cela dit, cela a le mérite de mettre les points sur les i: ce sont des physiciens. Pas des économistes.

Je pense qu’il y a des milliers de raisons d’attaquer les traders (si ce sont bien les “spécialistes” visés). Mais qu’on ne les prenne pas pour plus bêtes qu’ils ne le sont (je le dis d’autant plus volontiers que j’en ai formé un paquet, à l’ENSAE ou à Polytechnique, au cours des 10 dernières années).

Avant de revenir sur la petite histoire de martingale financière (tel que le résume le titre de l’article) faisons un court détour par celle brillamment racontée dans Mansuy (en ligne sur http://math.harvard.edu/~ctm/… et traduit en anglais dans http://jehps.net/Mansuy.pdf) sur l’origine de la notion de martingale. D’un point de vue mathématique (disons comme propriété des processus stochastiques), il faut remonter aux travaux de Joseph Bernstein, Paul Lévy, Émile Borel, et surtout Joseph Leo Doob, au milieu au XXième siècle. Cela dit, le mot est plus ancien: il entre dans le dictionnaire de l’Académie Française en 1762. Jouer à la martingale, c’est jouer toujours tout ce que l’on a perdu.  Mais quelques années plus tot, on pouvait déjà trouver le mot sous la plume de l’Abbé Prévost (le jeu qu’il décrit en 1750 comme variante du jeu du pharaon est aussi appelé martingale d’Alembert) Cela dit, l’origine que je préfère est de relier le mot martingale à l’expression provençale “jouga a la martegalo” qui signifierait jouer de manière incompréhensible, absurde, comme l’évoque Frédéric Mistral dans Lou Tresor dòu Felibrige ou dictionnaire provençal-français. On retrouve une origine proche dans le dictionarie of the French and English tongues de Randle Cotgrave, datant de 1611, qui mentionne l’expression “à la martingale” avec le sens absurdly, foolishly, untowardly, grossely, rudely, in the homeliest manner. Il cite même l’usage de l’expression philosopher à la martingale (sans citer Bernard Henri Levy, mais il semble que ce soit l’idée). Cela dit, même dans Lapinot on parle de martingale,

Bref, à partir de ce comportement incompréhensible, voire absurde, des économistes vont définir une notion importante en finance, que l’on appellera efficience (qui est un mot dangereux car il évoque aux oreilles de tous les économiste une notion d’optimalité au sens de Pareto, on pourra relire Malkiel (2003) The Efficient Market Hypothesis and Its Critics sur ce point). L’idée d’utiliser une marche aléatoire pour modéliser le cours d’un actif est ancienne. Par exemple (je ne remontrais pas au XVIème siècle, mais http://e-m-h.org/history.html le suggère) on pourra relire les travaux de Bachelier datant du début du XXième siècle, mais comme le dit la légende, Bachelier a été peu lu, à l’époque (en tous les cas par des économistes). Publié quelques années plus tard, on pourra relire Cowles & Jones (1937), Some A Posteriori Probabilities in Stock Market ou surtout Samuelson (1965) Proof That Properly Anticipated Prices Fluctuate Randomly. Cet article est passionnant, si on prend le temps de le lire “There is no way of making an expected profit by extrapolating past changes in the future price, by chart or any other esoteric devices of magic or mathematics.” Paul Samuelson est d’une modestie remarquable pour un chercheur, “I have not here discussed what the basic probability distributions are supposed to come from. In whose mind are they ex ante? Is there any ex post validation of them? Are they supposed to belong to the market as a whole? And what does that mean? Are they supposed to belong to the “representative individual”, and who is he? Are there some defensible or necessitous compromises of divergent anticipations patterns? Do price quotations somehow produce a Pareto-optimal configuration of ex ante subjective probabilities? This paper has not attempted to pronounce on these interesting questions“. Il peut y avoir des bulles, des processus non Gaussien, un peu tout ce qu’on peut imaginer, le point important étant qu’il est impossible d’utiliser le passé pour prédire le futur, en finance.

Oui, depuis 50 ans, les économistes savent que la marche aléatoire peut être un bon modèle pour décrire le prix des actifs. Ou pour être plus précis (et rendre à Paul Samuelson ce que Paul Samuelson a dit le premier) les martingales. Mais il faudra surtout attendre les travaux d’Eugène Fama, en particulier Fama (1965) The Behavior of Stock Market Prices et Fama (1970) Efficient Capital Markets: A Review of Theory and Empirical Work, qui vont de manière définitive marquer le début de l’utilisation des martingales en finance. Eugène Fama y retient trois notions d’efficience (pour aller plus loin que le fameux “A market in which prices at any time “fully reflect” available information is called “efficient”“)

1. Expected Returns or “Fair Game” Models
2. Submartingale Models
3. Random Walk Model

(en 1965, l’efficience devait être reliée, pour Eugène Fama, à la notion de marche aléatoire, alors que Paul Samuelson utilisait déjà la notion de martingale). Ces modèles disent tous (on est en 1970) qu’il est inutile d’utiliser l’information passée ou présente “to predict the future in a way which makes expected profits greater than they would be under a naive buy-and-hold model“. N’en déplaise à Pierre Barthélémy , depuis presque 50 ans ce principe est énoncé dans les cours d’asset pricing. Bachelier disait la meme chose en 1900, “l’espérance mathématique du spéculateur est nulle” dans Théorie de la spéculation. Et cela sera confirmé quelques années plus tard par les résultats empiriques de Cowles et Jones. Un an plus tard sera publié Black (1971) Implications of the random walk hypothesis for portfolio management, et trois ans plus tard, Hagerman & Richmond (1973) Random Walks, Martingales and the OTC. La littérature économique va se multiplier pendant les années 70, à tel point qu’en 1978, Michael Jensen, alors professeur à Harvard écrivait “I believe there is no other proposition in economics which has more solid empirical evidence supporting it than the Efficient Market Hypothesis” (cité dans http://economist.com/14030296). Plus récemment, à la fin des années 80 était publié le remarquable LeRoy (1989) Efficient Capital Markets and Martingales, qui étoffait la critique initiée dans LeRoy (1976). Efficient capital markets: Comment. LeRoy et Samuelson ont été les premiers à parler de martingale pour modéliser les prix des actifs. On continue (avec les articles qui font référence, et que tous les étudiants qui ont suivi un cours d’asset pricing ont lu) ? L’année suivante était publié Lehman (1990) Fads, Martingales, and Market Efficiency. Les martingales sont des outils incroyablement riches. Ils ne sont pas équivalent à une marche aléatoire: on peut avoir une martingale, sans avoir de marche aléatoire (par exemple avec un processus ARCH, on pourra relire sur le sujet les notes de cours Predictability of Asset Returns, en particulier le premier paragraphe, ou en français, avec des cours de licence, par Francis Diener,  http://math.unice.fr/~diener/…). Pour ceux qui veulent moins de formalisme, une martingale (en simplifant outrageusement) c’est, comme l’explique Nicolas Poupon,

Mais ces articles posent essentiellement les bases de la théorie financière. Qu’en est-il des aspects empiriques ? Andrew Lo et Craig MacKinlay ont publié un livre, en 2001, qui recense plusieurs articles sur le sujet, A Non-Random Walk Down Wall Street. Qui remet en cause l’idée de marche aléatoire, moins celle de marginale. On pourra aussi penser à Beechey Gruen & Vickrey (2000) et leur étude The Efficient Markets Hypothesis: A Survey, en ligne sur la Reserve Bank of Australia. Enfin, Pour une méthodologie plus proche de celle évoquée dans l’article en ligne sur arxiv, on pourra lire Martingales, the Efficient Market Hypothesis, and Spurious Stylized Facts de Joseph McCauley, Kevin Bassler et Gemunu Gunaratne (si on souhaite aussi utiliser des propriétés de mémoire longue). Moralité ? Non, ce n’est pas nouveau que l’on teste l’hypothèse d’efficience et de martingale. Et si je cite des articles académiques, il faut  ajouter que les journalistes économiques connaissent également cette littérature. Il y a 10 ans, Justin Fox posait la même question que nos chercheurs dans un article dans Fortune intitulé Is The Market Rational ? No, say the experts. But neither are you–so don’t go thinking you can outsmart it. On pourra aussi relire Efficiency and beyond (mentionné auparavant) qui revient sur la notion d’efficience des marchés. On y retrouve que les spécialistes savent tout ca: “on such ideas, and on the complex mathematics that described them [i.e martingales], was founded the Wall Street profession of financial engineering“. Enfin, en 2009, Richard Thaler faisait quelques rappels sur la difficulté d’interpréter correctement l’hypothèse d’efficience des marchés dans Markets can be wrong and the price is not always right. Pour conclure sur l’histoire des martingales en finance, en 2010, Fama répondait d’ailleurs de manière très claire à la question dans une interview – interview with Eugene Fama – lorsque John Cassidy lui demandait “the fundamental insight of the efficient market hypothesis [is] that you can’t beat the market ?“, et qu’il répondait “Right—that’s the practical insight. No matter what research gets done, that one always looks good“.

Pour revenir sur l’article de Pierre Barthélémy, je trouve qu’écrire “le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre” est incroyablement méprisant envers des centaines, pour ne pas dire des milliers de chercheurs en mathématiques financières, et en économétrie de la finance. Mais rassure toi Pierre, je vais continuer à lire tes chroniques. C’est juste que ma grand mère est encore abonnée au Monde, et j’aimerais bien qu’elle continue de croire que je peux sauver l’humanité (et non pas que je refuse d’admettre ce que des chercheurs italiens viennent d’établir, même si cela confirmerait – à ses yeux – que je peux etre têtu comme une mule quand je veux).

Maintenant, pour conclure sur l’article mentionné, et sur l’analyse faite, on ne peut pas dire que les conclusions soient révolutionnaires. Disons, pour quiconque ayant lu un peu de littérature économique dans les 50 dernières années. En fait, pour être honnête, je pense qu’on peut dire que Alessandro Pluchino et Andrea Rapisarda se sont amusé (on inclura Biondo et Helbing). Sans plus de prétention. Pour ceux qui ne se souviennent pas, ils avaient eu un IGNobel en 2011, pour leur étude sur organizations would become more efficient if they promoted people at random. J’avais aussi beaucoup aimé leur article sur les aspects computationnels du principe de Peter, en ligne sur http://arxiv.org/0907.0455. Quand je lis leurs études, j’ai l’impression qu’ils s’amusent. Et le terme n’est absolument pas péjoratif dans ma bouche, loin de là. C’est un peu ce que disait Lucile Quillet la semaine passée dans http://etudiant.lefigaro.fr/…  lorsqu’en mentionnant une étude de Baptiste Coulmont, elle écrivait “c’est le jeu auquel c’est amusé Baptiste Coulmont sur son blog“. Je ne sais pas si le terme était voulu, si c’était basé sur le fait qu’un blog, c’est fait pour jouer… mais j’ai toujours revendiqué que je faisais de la recherche pour mon plaisir, pour m’amuser. Et comme me le rappelait Baptiste, le mot schola signifie école, en latin, et son étymologie grecque est le mot σχολή, signifiant loisir.

# Rationality, and MS Excel (and other calculators)

This morning, Mathieu had a nice experience in his course on computational method in actuarial science. But let us start with some mathematical formal definitions.

First, recall that $y^x$ is – somehow – a standard expression. No one should be surprised to see such an expression. Generally (as explained in http://en.wikipedia.org/… ), this function is defined only when $y\in\mathbb{R}_+$. The idea is that the definition of $y^x$ is that

$y^x = \exp\left(x\log[y]\right)$

And it is a definition. Such a function exists only if $y\in\mathbb{R}_+$ (maybe excluding $0$). This would be a standard definition in real-analysis.

Now, this ‘power’ function appears also in complex analysis, when dealing with unit roots. From instance, if  $z=y^{\frac{1}{k}}e^{i \frac{2n\pi}{k}}$, where $y\in\mathbb{R}_+$ and $k\in\mathbb{N}_\star$, for some $n\in\mathbb{N}$, then $z^k=y$. Thus, in complex-analysis it might be more complex to define properly $y^x$ since it might not be unique. But we can relate (sometimes, when $x$ is the inverse of an integer, or maybe a rational number ?) with roots of polynomial functions. So far, nothing new…

Let us get back to Mathieu’s problem. Actually, in his course, he wanted to compute $(-8)^{\frac{1}{3}}$. With a French version of Excel, entering

you do get $-2$. If you look at the ‘help’ window, you have some more details

It looks like this hat function can be used to define objects such as $y^x$. But with

you get

(meaning that this is a problem…). It is also possible to use the power (puissance in French) function of Excel,

Here, you also get

The weird part here is that, in the ‘help’ window, you can read that this power function can be used with any number in $\mathbb{R}$.

Another point… what about $(-8)^{\frac{2}{3}}$ ? Somehow, it is just the square of the previous one (with the fraction)… Here, typing

you get

(similarly with the power function). So clearly, it is not that simple to use this power function. Now, if you use Google (which is now my new online calculator when I am in class, when I cannot use R), if the power is a fraction (or to be more specific the inverse of an integer), then it works as Excel

you get

But if you type (which should be close, from a continuity property of the power function)

you get

and similarly

On Wolfram Mathworld, enter

Mathematica does recognize that we try to deal with unit roots: the result is here

with – as expected – a numerical approximation

With Matlab, Mathieu did obtain the same as Mathematica (its decimal approximation). And to conclude, with R, Mathieu did obtain

> (-8)^(1/3)
[1] NaN
> (-8)^(.333333333333333)
[1] NaN

So for R, you cannot use this hat function on negative numbers.

Now, how can we interpret those outputs ?

1) My understanding is that clearly, with MS Excel, $x^{ab}\neq \left(x^a\right)^b$since

$(-8)^{\frac{2}{3}}\neq \left((-8)^{\frac{1}{3}}\right)^2$

which is problematic. For instance, in insurance, with monthly discounts, we do have functions like $u^{\frac{k}{12}}$. What if

$u^{\frac{k}{12}}\neq \left(u^{\frac{1}{12}}\right)^k$

2) The problem comes – probably (MS Excel is not an open software, so it might be hard to check) –  from the fact that $y^{\frac{1}{n}}$ is interpreted as an inverse of a (possibly) bijective function. To be more specific, $x=y^{\frac{1}{n}}$ means that $x^n=y$. When $n$ is an odd integer, then (in real-analysis) there is a unique inverse, and thus, $y^{\frac{1}{n}}$ is uniquely defined, since $x\mapsto x^n$ is a bijective $\mathbb{R}\rightarrow\mathbb{R}$ function. This is what MS Excel (and Google) is doing: $x\mapsto x^3$ is a bijective $\mathbb{R}\rightarrow\mathbb{R}$ function, so $(-8)^{\frac{1}{3}}$ means that we need to find the unique (real) value $x$ such that $x^3=-8$. Thus, somehow, it makes sense to return $-2$.

3) There is still a problem with Google, and Mathematica. That is fine to return unit roots in $\mathbb{C}$. But how comes there is only one value ? I mean, yes $1+\sqrt{3} \ i$ is a possible answer, since

$(1+\sqrt{3} \ i)^3=-8$

but one can also observe that , and similarly, $(-2)^3=-8$ and

$(1-\sqrt{3} \ i)^3=-8$

One can check with

With R, since we do not deal with power function here, but with roots, if we want to find $x$ such that $x^3=-8$, the function is

> polyroot(c(8,0,0,1))
[1]  1+1.732051i -2+0.000000i  1-1.732051i

Which is different… Weird isn’t it ?

# From Simpson’s paradox to pies

Today, I wanted to publish a post on economics, and decision theory. And probability too… Those who do follow my blog should know that I am a big fan of Simpson’s paradox. I also love to mention it in my
econometric classes. It does raise important questions, that I do relate to multicolinearity, and interepretations of regression models, with multiple (negatively correlated) explanatory variables. This paradox has amazing pedogological virtues. I did mention it several times on this blog (I should probably mention that I discovered this paradox via Marco Scarsini, who did learn me a lot of things, in decision theory and in probability). For those who do not know this paradox, here is an example that Marco gave in one of his talk, a few years ago. Consider the following statistics, when healthy people entered in some hospital

 hospital total survivors deaths survival rate hospital A 600 590 10 98% hospital B 900 870 30 97%

while, when sick people entered in the same hospitals

 hospital total survivors deaths survival rate hospital A 400 210 190 53% hospital B 100 30 70 30%

Somehow, whatever your health situation, you should choose hospital A. Now, if we agregate

 hospital total survivors deaths survival rate hospital A 1000 800 200 80% hospital B 1000 900 100 90%

i.e. without any doubts, people should choose hospital B.

and

while

With symbolic notations, one can have at the same time

and

with also

as shown on the graph below

There should be connection between Simpson’s paradox and the ecological fallacy (which is an issue I recently discovered and that I found extremely interesting, related again to difficulties of interpreting
regressions). But that’s another story. My point today is that Colin Blyth did mention another nice paradox, that is related, this time, to stochastic orderings. The idea is the following. Consider the three spinners drawn below (imagine some arrows in those circles)

• spinner A: no matter where the arrow stops, the gain is 3,
• spinner B: 56% chances to gain 2, 22% chances to gain 4, and 22% chances to gain 6,
• spinner C: 51% chances to gain 1, 49% chances to gain 5.

Instead of spinners, it is also possible to consider three different lotteries,

You play against a friend, you pick a spinner, while the friend picks another. Everyone flick his arrow, the highest number wins (no matter the difference). Let us compute the odds. First case, A against B, from
A’s perspective

 B-2 B-4 B-6 A-3 56% +1 win 22% -1 lose 22% -3 lose

In that case, A has 56% chance of beating B. Second case, A against C, from A’s perspective,

 C-1 C-5 A-3 51% +1 win 49% -2 lose
In that case, A has 51% chance of beating C. Third (an final) case, B against C, from B’s perspective. Assuming independence between the spinners, joint probabilities can easily be computed,
 C-1 C-5 B-2 28.56% +1 win 27.44% -3 lose B-4 11.22% +3 win 10.78% -1 lose B-6 11.22% +5 win 10.78% +1 win
In that case, B has 61.78% chance of beating C. So, if we try to summarize,
• A is the best choice, since it beats both with – always – more than 50% chance,
• C is the worst choice, since it is beaten by both with – always – more than 50% chance,
Now, assume that you play not against one friend, but two friends. An everyone picks a different spinner. Let
us compute the odds, one more time. First case, A against B and C, from A’s perspective
 B-2 C-1 B-2 C-5 B-4 C-1 B-4 C-5 B-6 C-1 B-6 C-5 A-3 28.56% +1 win 27.44% -2 lose 11.22% -1 lose 10.78% -1 lose 11.22% -3 lose 10.78% -3 lose
In that case, A has 28.56% chance of beating B and C. Second case, B against A and C, from B’s perspective,
 A-3 C-1 A-3 C-5 B-2 28.56% -1 lose 27.44% -2 lose B-4 11.22% +1 win 10.78% -1 lose B-6 11.22% +3 win 10.78% +1 win
In that case, B has 33.22% chance of beating A and B.Third (an final) case, C against A, from C’s perspective,
 A-3 B-2 A-3 B-4 A-3 B-6 C-1 28.56% -2 lose 11.22% -3 lose 11.22% -5 lose C-5 27.44% +2 win 10.78% +1 win 10.78% -1 lose

In that case, C has 38.22% chance of beating A and B. So, if we try to summarize, this time

• C is the best choice, since has (strictly) more than 1/3 chances to win, which the highest probability
• A is the worst choice, since has (strictly) less than 1/3 chances to win, which the lowest probability

Odd isn’t it ? Now, is there an interpretation of that paradox ? Yes, Martin Gardner, in his paper on induction and probability, mentioned the case of drug testing. The value we had with the spinner is the health level, rated from 1 to 6. Thus, taking drug A, you always get an average health level of 3. With drug C, on the other hand, you get either very sick (level 1) or very well (level 5). Consider now a doctor who wants to maximize the patient’s chance of being well. If only pills A and C are available, then the doctor should choose A. This is what we’ve seen in the first part. Assume that now a company delivers a third pill, called drug B. Then the doctor should find C more interesting…. Odd, isn’t it ?

Colin Blyth gave a more amusing application. Assume that you like to go to the restaurant, and you like get a dessert there. Dessert A – the apple pie – is the average one, with a standard level, that you rank 3 (on a scale from 1 to 6). Dessert C – the cheese cake – can either be awfull (ranked 1) or delicious (ranked 5). You’d better go for the apple pie if you want to maximize the probability of not being disappointed (i.e. maximizing your “best chance” according to Colin Blyth, but I guess it can be interpreted as regret minimization too). Now assume that dessert B – the blueberry pie – is available (with ranks given by the spinner). Then you should go for the cheese cake. I let you imagine the discussion that you can have, then, with your favorite waitress

– Hi Mr Freakonometrics, do you want a piece of apple pie ? (yes, actually she also comes frequently on my blog, and knows me from my pseudo…)

– Probably. But actually, I was wondering if you did have your blueberry pie today ?

– Yes, in fact we do….

– Great, in that case, I’ll go for the cheese cake.

She’ll probably think that I am freak… so I hope she’ll come and read my post, to understand that, actually, it does make a lot of sense to go for what was supposed to be my worst case.

# Sorting rows and colums in a matrix (with some music, and some magic)

This morning, I was working on some paper on inequality measures, and for computational reasons, I had to sort elements in a matrix. To make it simple, I had a rectangular matrix, like the one below,

> set.seed(1)
> u=sample(1:(nc*nl))
> (M1=matrix(u,nl,nc))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    7    5   11   23    6   17
[2,]    9   18    1   21   24   15
[3,]   13   19    3    8   22    2
[4,]   20   12   14   16    4   10

I had to sort elements in this matrix, by row.

> (M2=t(apply(M1,1,sort)))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    5    6    7   11   17   23
[2,]    1    9   15   18   21   24
[3,]    2    3    8   13   19   22
[4,]    4   10   12   14   16   20

Nice, elements are sorted by row. But for symmetric reasons, I also wanted to sort them by column. So from this sorted matrix, I decided to sort elements by column,

> (M3=apply(M2,2,sort))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    3    7   11   16   20
[2,]    2    6    8   13   17   22
[3,]    4    9   12   14   19   23
[4,]    5   10   15   18   21   24

Nice, elements are sorted by column now. Wait… elements are also sorted by row. How comes ? Is it some coincidence ? Actually, no, you can try…

> library(scatterplot3d)
> nc=6; nl=5
> set.seed(1)
> u=sample(1:(nc*nl))
> (M1=matrix(u,nl,nc))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    8   23    5   30   10   15
[2,]   11   27    4   29   21   28
[3,]   17   16   13   24   26   12
[4,]   25   14    7   20    1    3
[5,]    6    2   18    9   22   19
> M2=t(apply(M1,1,sort))
> M3=apply(M2,2,sort)
> M3
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    3    7   14   19   22
[2,]    2    6    9   15   20   25
[3,]    4    8   10   17   23   26
[4,]    5   11   16   18   24   29
[5,]   12   13   21   27   28   30

or use the  following function is two random matrices are not enough,

> doublesort=function(seed=2,nl=4,nc=6){
+ set.seed(seed)
+ u=sample(1:(nc*nl))
+ (M1=matrix(u,nl,nc))
+ (M2=t(apply(M1,1,sort)))
+ return(apply(M2,2,sort))
+ }

Please, feel free to play with this function. Because this will always be the case. Of course, this is not a new result. Actually, it is mentioned in More Mathematical Morsels by Ross Honsberger, related to some story on marching band. The idea is simple: consider a marching band, a rectangular one. Here are my players

> library(scatterplot3d)
> scatterplot3d(rep(1:nl,nc),rep(1:nc,each=nl), as.vector(M1),
+ col.axis="blue",angle=40,
+ col.grid="lightblue", main="", xlab="", ylab="", zlab="",
+ pch=21, box=FALSE, cex.symbols=1,type="h",color="red",axis=FALSE)

Quite messy, isn’t it ? At least, this is what the leader of the band though, since some tall players were hiding shorter ones. So, he brought the shorter ones forward, and moved the taller ones in the back. But still on the same line,

> m=scatterplot3d(rep(1:nl,nc),rep(1:nc,each=nl), as.vector(M2),
> col.axis="blue",angle=40,
+ col.grid="lightblue", main="", xlab="", ylab="", zlab="",
+ pch=21, box=FALSE, cex.symbols=1,type="h",color="red",axis=FALSE)

From the leader’s perspective, everything was fine,

> M=M2
> for(i in 1:nl){
+ for(j in 1:(nc-1)){
+ pts=m$xyz.convert(x=c(i,i),y=c(j,j+1),z=c(M[i,j],M[i,j+1])) + segments(pts$x[1],pts$y[1],pts$x[2],pts$y[2]) + }} But someone in the public (on the right of this graph) did not have the same perspective. > for(j in 1:nc){ + for(i in 1:(nl-1)){ + pts=m$xyz.convert(x=c(i,i+1),y=c(j,j),z=c(M[i,j],M[i+1,j]))
+ segments(pts$x[1],pts$y[1],pts$x[2],pts$y[2])
+ }}

So the person in the audience ask – one more time – players to move, but this time, to match with his perspective. Since I consider someone on the right, some minor adjustments should be made here

> sortrev=function(x) sort(x,decreasing=TRUE)
> M3b=apply(M2,2,sortrev)

This time, it is much bettter,

> m=scatterplot3d(rep(1:nl,nc),rep(1:nc,each=nl), as.vector(M3b),
+ col.axis="blue",angle=40,
+ col.grid="lightblue", main="", xlab="", ylab="", zlab="",
+ pch=21, box=FALSE, cex.symbols=1,type="h",color="red",axis=FALSE)

And not only from the public’ perspective,

> M=M3b
> for(j in 1:nc){
+ for(i in 1:(nl-1)){
+ pts=m$xyz.convert(x=c(i,i+1),y=c(j,j),z=c(M[i,j],M[i+1,j])) + segments(pts$x[1],pts$y[1],pts$x[2],pts$y[2]) + }} but also for the leader of the marching band Nice, isn’t it ? And why is this property always valid ? Actually, it comes from the pigeonhole theorem (one more time), a nice explanation can be found in The Power of the Pigeonhole by Martin Gardner (a pdf version can also be found on http://www.ualberta.ca/~sgraves/..). As mentioned at the end of the paper, there is also an interpretation of that result that can be related to some magic trick, discussing – in picture – a few month ago on http://www.futilitycloset.com/… : deal cards into any rectangular array: Then put each row into numerical order: Now put each column into numerical order: That last step hasn’t disturbed the preceding one: rows are still in order. And that’s a direct result from pigeonhole theorem. That’s awesome, isn’t it ? # Museum experiences, in the City This Saturday, we had two interesting museum experiences, with the kids. Kids are 10, 7 (and a half as she keeps saying) and 2 (and a half, too). In the morning, we went to the MoMaths (which is the pun which stands for Museum of Mathematics, and sounds like MoMA) see e.g. nytimes.com/2012/12/14/arts/…, which opened mid-December. There were a lot of exhibitions, to illustrate all kinds of mathematical concepts (as described in http://businessinsider.com/…). From a design point of view, what I did prefer was undoubtly the packing problem table. You pick up -say – five disks, and you have to pack them in – say – a triangle (an equilateral one). When you move the disks, the computer locates them and instantaneouldy computed the area of the smallest triangle which contains all of them: you have your own area, as well as the smallest one (known, so far). I wanted to see if kids had more intuition than I did, but even if we know the smalled area, unfortunately, we cannot see how to get it, in the museum. Of course, online, one can easily find it, e.g. for circles for triangles, or for squares Nice and fun, isn’t it ? (pictures are from http://mathworld.wolfram.com/…) About people around, you could hear two kinds of talks in the museum: “look, it’s fun, when you run, the ball is following you” (kids) or “oh, nice, as long as people move, the algorithm removes automatically nodes in the connected graphs” (former kids). Even some people decided to use sheets of paper to prove that they were correct, And yes, everyone had fun, even my youngest daugther, who’s only 2 (and a half, I know). We spent two hours at the MoMaths, and the kids did not want to leave. I have to admit that I was a bit relieved to leave, finally: visiting the museum was like those social events in conferences. Actually, people there were clearly people that I keep seing at work, and in conferences. But there were kids, too, so it was much more convivial. People were also drawing, like in real museums… or like during talks, in conferences. It was work, with fun (and the family). Like my blog, somehow… Then, we went to the Guggenheim Museum. One has to admit that this is another experience. The building itself is amazing. But here, it was a more standard museum experience… Except perhaps while we where visiting the Gabriel Orozco exhibition (which was, actually, extremely interesting and reminded me Francis Alÿs ‘s work presented at MoMA last year, http://francisalys.com/…) There was a guard, in a corner, and when we got too close, he asked us to move back. Trying to understand what was going on, we saw a small stone, behind him (of course “no picture, please“). The guard was waiting for a museum conservator to assess if the stone was from a shoe of a visitor, or a part of the exhibition… The conservator came look clearly at the stone, as well as others, trying to choose (if possible with conviction) whether this was a piece of art (part of the exhibition), or just a vulgar stone. True story… We finally left after one hour (my daughter had finished her kids activity, i.e. drawing like Picasso, comparing Picasso and Velasquez’s Meninas, and counting nuances of grays). Now, let us compare those two experiences. For instance, the visitors. At MoMaths, it was half kids, half math professors (from a personal guess), while at Guggenheim, it was more like one third tourists (with their NY Yankees cap, just to pretend they aren’t tourists… which is exactly what I do), and two third retired upper-class New Yorkers (including some much youger models in the arms of golden boys). At MoMaths, most of them were English speakers (and they had to, beacuse it was quite difficult to get information in another language, unfortunately), but at Guggenheim, I heard a lot of French, Italian, even Dutch… The public in those two museums were completely different. When discussing with kids afterwards, we wanted to know which museum they did prefer. Without any doubts, the answer was MoMaths. And when we asked them in which museum the freakiest people were, they both said the Guggenheim museum. Maybe I am not really a freak, after all, at least in my kids eyes… # Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 2) Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia, > capital=read.table( + "http://freakonometrics.free.fr/US-capital.csv", + header=TRUE,sep=";") > capital=capital[-c(2,11),] On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales, > library(maps) > library(maptools) > library(sp) > library(TSP) > Lcap=paste(as.character(capital[,1]), + as.character(capital[,2]),sep=", ") > M=as.matrix(USCA312) > L=labels(M)[[1]] > Lcap=Lcap[which(Lcap%in%L)] > k=which(L%in%Lcap) > M=M[k,] > M=M[,k] > listeUSA=TSP(M,L[k]) Et cette fois on est prêt, > tour=solve_TSP(listeUSA, method = "nn") on a lancé l’algorithme, et on peut faire un dessin, > COORD=data.frame(USCA312_coords[k,]) > COORD=COORD[as.numeric(tour),] > map('state',fill=TRUE,col="yellow") > lines(COORD$coords.x1,COORD$coords.x2, + lwd=3,col="blue") > points(COORD$coords.x1,COORD\$coords.x2,
+ pch = 19, cex = 1.2, col = "red")

pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet,

voire encore un autre si on lancer une nouvelle fois l’algorithme,

C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante

autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard,

Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité)

>  HAZARD=OPTIMAL=RATIO=rep(NA,500)
>  for(s in 1:500){
+   HAZARD[s]=tsp.longueur(M,sample(1:nrow(M)))
+   tour=solve_TSP(listeUSA, method = "nn")
+   OPTIMAL[s]=tsp.longueur(M,as.numeric(tour))
+   RATIO[s]=HAZARD[s]/
+            OPTIMAL[s]
+ }

Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant

> tour=solve_TSP(listeUSA, method = "nn",
+ control=list(start=41))

Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths…