Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Le prochain qui me parle du R2…

Depuis que j’ai repris la recherche, j’ai pris l’habitude de débrancher mon téléphone. Mais lundi matin, comme j’avais un coup de fil à passer, j’ai du le rebrancher. J’ai alors reçu un appel d’un élève qui m’a agacé (et auquel j’ai eu la lâcheté de raccrocher en expliquant que je faisais autre chose – ce qui était néanmoins partiellement vrai) mais qui est assez symptomatique, “bonjour, j’ai fait une régression et j’ai un R2 de 0,312 et je voulais savoir si c’était bien ?” (ou “si c’était beaucoup ?”, ou “si c’était assez ?”, ou un truc du genre car j’avoue que j’avais effectivement la tête ailleurs).
Pourtant je croyais que mon cours avait précisément réussi à éviter de parler du R2. Je vais donc reprendre ici quelques points. J’espérais trouver un site internet de l’association des statisticiens qui en ont marre du R2 (ou d’opposants à la sectes des adorateurs du R2, car si j’ignore qui en est le gourou, je trouve que cette secte regroupe beaucoup d’adorateurs) malheureusement je n’ai rien trouvé de concluant sur internet. Donc je vais m’y mettre… même je risque de faire un très long billet,

  • définition et interprétation du R2


Un “bon” modèle permettra d’obtenir des estimations “proches” des valeurs observées. Sur la représentation dans l’espace des variables (comme le montre la figure ci-dessus) la qualité peut être évaluée par l’angle θ. Cet angle est compris entre -90˚et 90˚. Un angle proche de -90˚ou de 90˚indique un modèle de mauvaise qualité. Alors qu’un angle proche de 0˚ correspond à un bon modèle. Le cosinus carré de θ est donc une mesure possible de la qualité du modèle et cette mesure varie entre 0 et 1.
En utilisant le théorème de Pythagore nous permet d’écrireou encoreoù SCT (respectivement SCE et SCR) désigne la somme des carrés totale (respectivement expliquée par le modèle et résiduelle). Le coefficient de détermination R2 est défini parc’est-à-dire que le R2 est souvent interprété comme le “pourcentage de la variabilité de la variable endogène expliquée par le modèle”.
Si on reprend une page du livre de Pierre-André Cornillon et Éric Matzner-Løber, on a les définitions suivantes

  • que mesure le R2 ? quelques exemples…

Comme je le dis auparavant, le R2 nous dit si la courbe de régression (pour prendre un modèle plus général que le cas linéaire) est proche des observations, mais ne nous dit nullement si le modèle est “bon“.
Pour faire simple, avec des observations individuelles, le R2 sera souvent assez “faible“, alors qu’avec des séries temporelles (en particulier des séries intégrées), le R2 sera souvent “plus élevé“. Dans ce dernier cas, l’explication heuristique est qu’expliquer une série croissante par une autre série croissante marche souvent “bien“.Sur l’exemple ci-dessus, on obtient – en ajustant un modèle linéaire – respectivement
Multiple R-squared: 0.3372,
Multiple R-squared: 0.6414,
Sur l’exemple ci-dessus, on obtient – en ajustant un modèle polynôme de degré 4 – respectivement
Multiple R-squared: 0.4968,
Multiple R-squared: 0.987,

Sur l’exemple 3, l’ajustement est “bon” au sens où il serait difficile de faire mieux, c’est simplement qu’il est très bruité. Avec un R2 de 0,5, l’exemple 3 me paraît proposer un “meilleur” ajustement que l’exemple 2, alors que le R2 est de 0,65. Dans le premier cas, l’ajustement semble bon, mais très bruité, alors que pour le second, le modèle est plutôt mauvais, mais le bruit est relativement faible.

  • le R2 n’est qu’un coefficient de corrélation, ou presque (et donc possède tous les défauts d’un coefficient de corrélation)

Si on reprend le premier point, il est possible d’écrire le R2 (dans un modèle avec une unique variable explicative) comme le carré du coefficient de corrélation (au sens de Pearson) entre la variable endogène et la variable explicative (plus généralement, ça serait entre les observations et la prédiction). Regardons par exemple sur la base de données cars,
> data(cars)
> cor(cars)
speed dist
speed 1.0000000 0.8068949
dist 0.8068949 1.0000000

Si on regarde la valeur du R2 on peut écrire
> summary(lm(dist~speed,data=cars))$r.squared
[1] 0.6510794
> sqrt(summary(lm(dist~speed,data=cars))$r.squared)
[1] 0.8068949

J’avais dit dans un précédant billet (ici) tout le mal que je pense du coefficient de corrélation, qui n’est pas une mesure de dépendance (ou une mesure de concordance au sens où l’avait défini Marco Scarsini, ici). Pour cette raison, je trouve qu’un R2 qui vaut 0,312 nous apprend aussi peu de chose qu’un coefficient de corrélation de 0,217 ou 0,013.

  • mais le R2 n’est pas exactement le carré d’une quantité (et donc peut être négatif, par exemple)

Formellement, le R2 est le carré du cosinus d’un angle, si on reprend l’interprétation géométrique de la régression linéaire. Mais si on regarde la version empirique, on peut être (désagréablement) surpris… En particulier avec Excel. Ce point est discuté sur de nombreux forum (ici ou ), je ne vais donc pas en rajouter.

  • le R2 augmente en rajoutant des variables explicatives (ce qui le rend difficilement utilisable pour comparer des modèles)

Ce point est souvent connu par tout le monde, c’est pour cela qu’on présente toujours le R2 ajusté après avoir présenté le R2.

  • le R2 ne peut pas être utiliser pour comparer des modèles pour lesquels la variable endogène change

Mais ceci n’est pas qu’un problème avec le R2, car il est très difficile de trouver un critère de choix suffisamment robuste. J’avais déjà abordé ce point lors de mon cours d’assurance non-vie à l’ensae il y a quelques années. J’avais vu des élèves essayer de comparer deux modèles (pour tarifer), à savoir

  • une régression gamma, i.e.
  • une “régression lognormale”, c’est à dire que

le soucis est qu’il est difficile de comparer des deux modèles. En particulier, tous les élèves qui se posaient cette question à l’aide d’un critère de type AIC concluaient toujours à la supériorité de la loi lognormale, ou pour être un peu pointilleux, sur un modèle gaussien après une transformation logarithmique.
Dans le cas des modèles Gaussiens, c’est aussi la conclusion à laquelle arrive Valérie Mignon dans son livre d’économétrie, à savoir que si l’on cherche à comparer

  • un modèle additif, i.e.
  • un modèle multiplicatif, i.e.

on retient “toujours” le modèle multiplicatif si on retient le modèle qui a le meilleur R2. Je n’ai pas de raison a priori d’aller dans ce sens, mais je retiens que juger de la supériorité d’un modèle à l’aide uniquement de ce critère me paraît vraiment douteux.
L’exemple ci-dessous représente les deux régressions, avec un modèle linéaire versus un modèle exponentiel (i.e. avec comme variable endogène le logarithme).
Dans ce cas précis,
> REG.LINEAIRE=lm(Y~X)
> REG.MULTIPLICATIF=lm(log(Y)~X)
> summary(REG.LINEAIRE)$r.squared
[1] 0.8688
> summary(REG.MULTIPLICATIF)$r.squared
[1] 0.8654
Afin de répondre à cet exemple précis, je renvoie au test de Box-Cox (que j’avais abordé ici) qui vise précisément à comparer ces deux modèles. Ici, on conclue à la supériorité d’un modèle multiplicatif, mais ça n’a rien à voir avec la valeur du R2.

  • le R2 peut être changé artificellement (en réécrivant le modèle)

Je reprendrais un exemple tiré d’une épreuve du concours interne d’administrateur de l’insee (qui doit avoir une vingtaine d’années). Considérons un modèle assez classique en économie, où on essaye de lier taux d’intérêt et taux d’inflation. Je renvois à mes slides du cours d’économétrie où l’exemple est présenté, et mis en œuvre, ici (slides 39 à 42).

  • le R2 peut être changé artificellement (en supprimant les points aberrants)


Ce point est présenté dans certains livres, où il est expliqué qu’en enlevant des points (souvent les premières années quand on travaille sur des séries macroéconomiques) on peut améliorer le R2. Effectivement, si on enlève les quatre points les plus éloignés de la prédiction, on change le modèle, et le R2 est alors amélioré,
> summary(REG.TOTAL)$r.squared
[1] 0.6510794
> summary(REG.SUBSET)$r.squared
[1] 0.7539631

Je croyais là aussi que mon cours essayait d’expliquer l’importance (et le traitement) de ces points aberrants.
Pour conclure, je reprendrais ce que j’essaye d’expliquer dans mon cours d’économétrie (ici ou là), à savoir qu’au lieu de se contenter de regarder un R2, il vaut mieux passer un peu de temps à regarder quelques graphiques de diagnostique. Et pour comparer des modèles, je préfère toujours les critères AIC ou BIC au R2 (ou même au R2 ajusté). Quant à la normalité des résidus… je pense que je reviendrais là dessus aussi dans un billet les jours à venir….

Who wants to be a millionnaire ?

Le vrai moyen de gagner beaucoup est de ne vouloir jamais trop gagner et de savoir perdre à  propos.” (Fénélon – Les aventures de Télémaque)

Bon, pour ma part, j’avais tenté une autre solution vendredi dernier: tenter de gagner les 126 millions d’euros de l’Euromillion. Et comme je suis toujours en train de poster des billets sur mon blog, on peut en déduire que je ne suis pas l’heureux vainqueur (car oui, il parait qu’il y en a un… un très bon ami ami à moi… enfin… s’il le souhaite!).

Quelle ne fut pas me surprise de voir des collègues (on peut dire copains) matheux indignés par mon comportement, et me reprocher mon irrationalité*

Sur l’irrationalité des joueurs, je renvoie à de nombreux papiers de la littérature économique, en particulier le papier de Maurice Allais (ici), précisément intitulé Le comportement de l’homme rationnel devant le risque. Je dois avouer que, selon moi, la critique de Maurice Allais est très proche de ce que notent Daniel Kahneman et Amos Tversky (ici) sur la distinction fondamentale (et non symétrique) entre gains et pertes.

Bref, toute cette discussion est un prétexte pour parler d’un très joli papier que m’a fait découvrir Charles Figuières lors de mon exposé à Montpellier voilà deux semaines, de Graciela Chichilnisky sur Topology of Fear. Un preprint peut se télécharger ici, sa page internet (ici) regorge d’informations passionnantes, mais je ne sais toujours pas où ce brillant papier sera publié…. à suivre donc.

* pour aller un peu plus loin, je pourrais avouer ici que ma stratégie a été de jouer “1-2-3-4-5” puis “6-7”  (j’ai découvert que ce type de loto était un peu différent de celui que j’utilise dans mes exemples de probas, avec des étoiles). Tout d’abord, cette stratégie est aussi probable que n’importe quelle autre. Mais pour information, lors du 3ème tirage du Loto en France (en 1976), les numéros “1-2-3” sont sortis. Historiquement (voir ici), sur plus de 4850 tirages, on a eu deux numéros consécutifs (au moins) dans 2/3 des tirages (63%), et trois consécutifs dans 8% des tirages. On a vu 39 fois (0,8% des tirages), quatre numéros consécutifs. Mais jamais 5 numéros consécutifs, en presque 35 ans. Je laisse les plus motivés faire les calculs théoriques. Ma stratégie devrait être payante (si je jouais plus d’une fois tous les 5 ans) car peu de monde joue des numéros consécutifs. Si on compte le nombre de gagnants du 1er rang, il y a – en moyenne – 2,90 gagnants les tirages comprenant deux numéros consécutifs et sinon 4,03 gagnants. Si on va plus loin, il y a 2,48 gagnants les tirages comprenant trois numéros consécutifs et 3,39 sinon. Par contre, si on regarde avec 4 nombres consécutifs, un tombe sur un point aberrant: d’ordinaire, on a 3 gagnants du 1er rang, mais en septembre 1984, avec comme tirage “26-27-28-29-35-36-37” (4 et 3 consécutifs), il y a eu 64 gagnants du 1er rang, ce qui est presque le record…
D’ailleurs si on regarde les records (promis après j’arrête), le record fut en juin 1993, avec 103 vainqueurs, où 31 fut le plus grand chiffre, et le 13 sortit. Ensuite, en 1979, on a eu 96 vainqueurs avec des numéros pairs consécutifs, à savoir le tirage “2-8-10-12-14-16-18” (soit 6 chiffres pairs successifs). Ensuite arrive le tirage que je mentionnais… En fait, tous les tirages qui ont eu des très gros nombres de vainqueurs du 1er rang présentaient une régularité impressionnante. On retrouvera ainsi le tirage “4-5-7-14-15-17-20” par exemple, avec 63 vainqueurs… Bref, ma stratégie n’est peut être pas géniale finalement….

Cube COSO

En préparant mes slides pour la conférence ERM – enterprise risk management – pour la Banque de France (ici), je n’ai pu m’empêcher de chercher des références sur la formalisation du COSO (le fameux référentiel de contrôle interne du Committee Of Sponsoring Organizations of the Treadway Commission utilisé dans le cadre de la mise en place des dispositions relevant de la loi Sarbanes-Oxley aux Etats-Unis). On y trouve en particulier une description assez fournie des risques, avec

  •  l’identification des événements potientiels (tendances, événements passés)
  •  l’évaluation des risques (risque inhérent, risque résiduel) ,
  •  les réponses aux risques (catégorisation des types de réponses).

(que l’on retrouvait déjà évoqué dans un vieux billet, ici). De manière visuelle, on retrouve d’ailleurs souvent le fameux cube,

parfois plus connu dans sa version couleur,

Mais en cherchant un peu, j’ai retrouvé cette représentation dans un vieux papier de Matthias Haller dans les Geneva Papers on Risk and Insurance de 1976, dans un papier intitulé “les objectifs du risk management“,

Différencier (indéfiniment) les séries temporelles ?

Dans un mail, un étudiant qui finissait son projet de séries temporelles m’a posé une question simple et intéressante  (que je me permets de reprendre ici): “quand on cherche à stationnariser une série, on a souvent  besoin de différencier deux fois, et ça marche tout le temps“.
Effectivement, on peut toujours différencier une série, on finira bien par tomber sur une série stationnaire. Je retraduirais cette question sous la forme suivante “pourquoi cherche-t-on toujours à stationnariser les séries ?” ou encore “est-ce gênant de différencier une série ?“.
Pour la première question, la réponse est simple: les seules séries que l’on sache modéliser sont les séries stationnaires. Les autocorrélations se calculent sur des séries stationnaires, par exemple, et la non-stationnarité n’est pas une notion simple à définir (tout comme la non-indépendance entre variables aléatoires)…..etc.
Pour la seconde question, la réponse est simple: si on différencie une série pour mieux la modéliser, si on souhaite par la suite faire de la prévision, il conviendra d’intégrer la série modélisée. Or intégrer une série fait exploser l’intervalle de confiance. Autant faire ça sur un exemple simple, avec la série suivante, que l’on commence par supposer stationnaire, et que l’on prédit sur une trentaine de valeurs, avec un intervalle de confiance.
XXX

Mais si l’on suppose la série intégrée à l’ordre 1, on différencie puis on modélise par un processus stationnaire (c’est l’idée qui sous-tend les processus ARIMA, à savoir un processus ARMA intégré: en différenciant la série, on devrait tomber sur un processus ARMA. On peut rattacher ça à la notion de racine unité). La prédiction (avec l’intervalle de confiance) est alors présenté sur la série initiale. Le fait d’intégrer les erreurs fait que la variance de la prédiction augmente avec l’horizon.

Et si l’on continue, et que l’on différencie de fois la série avant de la modéliser, et que l’on intègre deux fois les erreurs, l’intervalle de confiance devient immense !
XXX

Bref, il convient de ne pas différencier si ce n’est pas vraiment indispensable ! L’outil pour vérifier que l’on n’a pas trop différencier est la fonction d’autocovariance inverse (c’est à dire la fonction d’autocovariance d’une série dont la densité spectrale serait l’inverse de la densité spectrale initiale*). Cette fonction fait partie des fonctions de base sous SAS, mais pas sous R… J’ai cherché un peu, mais sans succès. L’argument dans les forums est que cette fonction est redondante avec la fonction d’autocorrélation partielle. Et en effet, si la finalité est simplement de détecter les ordres AR ou MA d’une série stationnaire, alors effectivement, les deux fonctions s’utilisent de manière équivalente. Mais la fonction d’autocovariance inverse apporte plus d’information, en particulier afin de détecter si l’on n’a pas surdifférencié la série.
La sortie suivante présente la fonction d’autocorrélation (à gauche) et la fonction d’autocorrélation inverse (à droite) sur une série.
XXXX

On note que la série semble intégrée (en pratique, des autocorrélations très fortes et persistantes très longtemps se traduit par une suspission d’intégration). Si l’on différencie la série, on obtient les autocorrélogrammes suivants,
XXX

Et si l’on différencie une nouvelle fois, on obtient les graphiques suivants,
XXXX

Bref, l’autocorrélogramme de droite semble caractéristique d’une série intégrée: on dira alors que l’on a surdifférencié le modèle (“en intégrant la série, elle est toujours stationnaires”…). Je parle un peu de tout ça dans mon polycopié de Dauphine, tome 1.

* J’avais parlé dans un ancien billet (ici) de l’importance de la densité spectrale en série temporelles. L’idée est que si f est une densité spectrale, 1/f peut également l’être. Le lien entre la fonction d’autocorrélation et la densité spectrale est donné par des théorèmes de Kolmogorov ou Wiener. Pour plus de détails, il y a des compléments dans le chapitre 1.2 du polycopié que j’avais fait à Dauphine (ici). On peut d’ailleurs noter que la fonction d’autocorrélation partielle d’un processus ARMA(p,q) est la fonction d’autocorrélation d’un processus ARMA(q,p) obtenu en permutant les deux polynômes d’opérateurs retard.

Séminaire Probabilité et Statistique, UBO, Brest

Talk at the statistical seminar at the Université de Bretagne Occidentale, in Brest, Wednesday May 6th Tuesday May 5th, 14h (in  10 days), on “multivariate extremes. Slides can be found here.

The talk will give a detailed introduction on multivariate extremes and related concepts. Then the case of Archimedean copula will be fully described (following the paper with Johan Segers).

[04/05/2009]: some applications in risk management will be shown at the end of talk, as well as some news things on spatial correlation.

and in order to illustrate tail convergence of Archimedean copulas, I have uploaded two animations, with tail independence below,

with tail dependence (or asymptotic dependence),

Expliquer l’assurance vie et la retraite

En fouillant un peu dans mes cartons, j’ai enfin remis la main sur un article qui avait trôné dans mon bureau à l’ENSAE pendant plusieurs mois. De mémoire, je croyais qu’il s’agissait d’un article du Monde, mais après quelques recherches, je me suis rendu compte qu’il s’agissait d’un article paru fin 2003 dans le Canard Enchaîné, hebdomadaire pour lequel j’ai une immense admiration pour les informations (et surtout les jeux de mots), mais pour une fois, je serai assez critique.

En particulier, le petit encadré en bas à gauche m’avait interpellé,

Cet article parle du calcul des retraites – sujet ô combien important – et s’étonne, ou s’insurge, que le calcul d’une retraite soit “insensé” (c’est à dire selon mon dictionnaire “qui a perdu le sens“). Or précisément cette formule a du sens (tout du moins si l’on remplace x – qui est l’inconnue révélatrice de l’incompréhension face aux formules mathématiques dans l’imaginaire collectif – par un signe de multiplication). Les rapports des L sont les probabilités d’être en vie – et donc de pouvoir toucher effectivement une retraire, et on retrouve les facteurs d’actualisation devant. Bref, on retrouve une formule qui s’interprète “simplement” comme valeur actuelle probable des flux futurs. Le document initial peut se trouver ici (merci à Benoit Bellone pour la référence), il s’agit de la section 8.3 (“Garantir le niveau des retraites: Les conséquences financières pour les futurs retraités”) du Rapport économique, social et financier, Tome I  Perspectives économiques 2003-2004 et évolution des finances publiques. Une partie explique en détail la “neutralité financière d’un régime de retraite“, et ce sont ces formules qui semblent avoir inspiré l’article du Canard,

Je pense que cet article (pourtant paru dans un journal dont le lectorat est majoritairement de la classe supérieure élevée) est assez symptomatique du manque de culture scientifique des journalistes. Sans vouloir faire offense à l’esprit anticlérical du Canard, cet article m’a rappelé une réflexion du philosophe Calvin (ami de Hobbes, le tigre)

"Générer" du hasard, partie 5

La dernière méthode que  je présenterais est est celle des générateurs congruentiels linéaires. On retrouve l’idée de linéarité que l’on avait pas exemple dans le premier générateur proposé (ici), la notion de congruence est liée à des problèmes d’arithmétique et d’écriture de nombres dans les bases. On retiendra une fonction de récurrence de la forme suivante

X_{n+1} = (a cdot X_n + c) mod m

Pour que l’algorithme marche bien, on essaye de partir d’une valeur initiale correspondant à un nombre premier.
Prenons un premier exemple simple (inspiré par wikipedia). Si on prend des valeurs un peu au hasard, par exemple on obtient

  • avec X0 = 10, la suite : 10, 10, 10, 10, 10, …
  • avec X0 = 11, la suite : 11, 35, 123, 19, 235, 3, 91, 243, 203, 227, 59, 211, 171, 195, 27, 179, 139, 163, 251, 147, 107, 131, …
  • avec X0 = 12, la suite : 12, 60, 236, 28, 204, 252, 172, 220, 140, 188, 108, 156, 76, 124, 44, 92, 12, 60, 236, 28, 204, 252, 172, …

Bref, on notera déjà qu’il est effectivement intéressant de partir d’un nombre premier, car on boucle alors beaucoup moins vite. Ah oui, et pour se ramener sur l’intervalle [0,1], on divise par m. On a alors intérêt à avoir m assez grand.
Prenons un second exemple simple (inspiré du générateur proposé par Robert Sedgewick). On prend , et on part de 0.

[1]         0        1 31415822 40519862 31536702 84962343
[7]  59428603 54128063 38284723  4802584 19281465 53057765
[13] 47900065 67928365 55662665 20022965 84329266 25717386
[19] 95163907 67972648 42464008 74270568 69856328 96165288
[25] 74221448 24728808  5671369 13328949 75902129 98182909

La série a l’allure suivante,
xxxxEt si on regarde l’histogramme et la fonction d’autocorrélation sur les 10 000 premières valeurs, on a gardé l’uniformité (ou tout du moins on ne doit pas en être trop loin), et on n’a plus d’autocorrélation.

Parmi les générateurs ultra-classiques, mais mauvais, on retiendra le RANDU, datant des années 60, implémenté alors sur les IBM, X_{n+1} equiv (65539 times X_n) mod 2^{31} Mais cet algorithme présente de très gros défauts numériques. Je n’aborderais pas ici les problèmes de ce générateur (je l’avais fait dans le chapitre sur les simulations dans le tome 2 du livre avec Michel Denuit).
Bref, voilà ce qu’on peut raconter rapidement sur les méthodes de construction de nombres aléatoires, ou comme on dit généralement pseudo-aléatoires. Pour quelques réflexions sur l’utilisation de méthodes déterministes pour générer du hasard, je renvoie au très beau petit livre d’Ivar Ekeland intitulé Chaos. Et histoire d’insister un peu lourdement, ces billets (comme la plupart de ceux présent sur ce blog) n’ont pas grande prétention, si ce n’est présenter des concepts parfois compliqués “avec les mains” (et des dessins). Pour aller plus loin, on trouvera de très bons résumés de ces techniques sur le net, ici ou . Je laisse les plus curieux continuer. Sinon le chapitre 7 du Monte Carlo de Fishman propose une analyse poussée de ces méthodes. Et je voulais aussi noter que ce que j’ai proposé dans ces billets, c’est simplement de tirer des nombres aléatoires uniformément répartis sur [0,1]. Pour générer autre chose que de l’uniformité, je renvoie à la bible de Luc Devroye, non-uniform random variate generation, épuisé chez Springer, mais en ligne ici.

"Générer" du hasard, partie 4

Bon l’idée de récurrence les dernières fois (ici et ) n’était pas forcément mauvaise, mais telle qu’elle était présentée, ça ne marchait pas très bien… Une autre méthode peut être celle proposée par John von Neumann, parfois appelée middle square. On prend un nombre à 4 chiffres, on l’élève au carré, et on prend les 4 chiffres du “milieu“.  Par exemple à partir de 1234, dont le carré est 1522756, puis on prend le carré de 2275, c’est à dire 5175625… etc. On prend alors comme nombre cette suite de nombres divisés par 10000 (pour revenir entre 0 et 1). Sur cette exemple on note que l’on peut assez rapidement boucler, et retomber sur la même suite (ceci reste le problème central de tous les algorithmes récursifs).
XXX

Bref, on n’y est toujours pas… Pourtant, de loi, la suite des nombres peut  sembler aléatoires (en oubliant le fait que l’on s’est restreint à 4 chiffres),

[1]  0.2756 0.5536 0.7296 0.1616 0.1456 0.9936 0.4096 0.7216
[9]  0.7065 0.4225 0.5062 0.3844 0.6336 0.4896 0.7081 0.4056
[17] 0.1136 0.0496 0.6016 0.2256 0.9536 0.5296 0.7616 0.3456
[25] 0.3936 0.2096 0.3216 0.2656 0.4336 0.8008 0.8064 0.8096
[33] 0.5216 0.6656 0.2336 0.6896 0.4816 0.3856 0.8736 0.7696

Mais si on regarde cet exemple de plus prêt, on note que le chiffre 6 apparaît très fréquemment à la fin, par exemple…. Même si cette méthode est souvent présentée (par exemple dans la page de wikipedia consacrée aux méthodes de simulations, ici), elle marche beaucoup moins bien que les méthodes que j’avais présentées auparavant (ici et ) et que l’on cite assez peu.

"Générer" du hasard, partie 3

Pour poursuivre l’idée développée dans le billet précédant (ici), prenons une autre fonction de récurrence.

Toujours dans ces notes de cours (donné il y a un peu plus de 10 ans par Christian Gouriéroux), il y avait cette tent function que l’on retrouve sur le graphique de gauche. Si je fais tourner l’algorithme à partir de π-3, on obtient, avec une simple fonction R

[1]  0.14159265 0.28318531 0.56637061 0.86725877 0.26548246 0.53096491
[7]  0.93807017 0.12385966 0.24771932 0.49543864 0.99087728 0.01824545
[13] 0.03649090 0.07298179 0.14596358 0.29192717 0.58385434 0.83229132
[19] 0.33541736 0.67083471 0.65833057 0.68333886 0.63332228 0.73335543
[25] 0.53328913 0.93342173 0.13315654 0.26631308 0.53262615 0.93474770

Bref, ça marche bien là aussi (à condition là aussi d’avoir une écriture très longue de π). On peut là aussi montrer simplement que l’on génère une suite de nombre uniformément distribués sur [0,1], car cette loi est invariante par la fonction de récurrence (tout comme celle proposée par Whittle). Mais on n’a toujours pas d’indépendance….

“Générer” du hasard, partie 2

Une première piste (qui n’est pas de moi mais qui a été proposé par Peter Whittle en 1963, et que j’ai retrouvé dans des notes prises lors d’un cours de Christian Gouriéroux) est une idée simple de récurrence. On prend un nombre u entre 0 et 1 (si possible irrationnel, sinon ça n’a pas grand intérêt), et on considère la suite. Par exemple si on prend comme valeur initiale π-3, on obtient, avec une simple fonction R

[1]  0.14159265 0.28318531 0.56637061 0.13274123 0.26548246 0.53096491
[7]  0.06192983 0.12385966 0.24771932 0.49543864 0.99087728 0.98175455
[13] 0.96350910 0.92701821 0.85403642 0.70807283 0.41614566 0.83229132
[19] 0.66458264 0.32916529 0.65833057 0.31666114 0.63332228 0.26664457

Bref, ça marche bien (à condition numériquement d’avoir une écriture avec un nombre infini de décimales).  Le plus simple est de passer en base 2 (on peut trouver π en base 2 sur le net, par exemple ici). Bref, on peut obtenir la série suivante,
Plus généralement sur le caractère aléatoire des décimales de π, je renvoie aux publications de David H. Bailey (ici).
Bref, cette méthode pourrait marcher car on génère effectivement des nombres uniformément sur l’intervalle [0,1]… sauf que les tirages ne sont pas du tout indépendants. On observe un beau processus autorégressif AR(1), avec une autocorrélation  en 2-h.
XXBref, en l’état cette méthode n’est pas vraiment géniale pour générer du “hasard“… à suivre donc.

“Générer” du hasard, partie 1

J’ai déjà mis plusieurs billets sur mon blog afin d’évoquer les méthodes de Monte Carlo (ici ou là), car de mon point de vue, ce sont des techniques fondamentales en économétrie par exemple, si on travaille sur des petits échantillons. Mais pour l’instant, j’avais toujours éludé le problème délicat du générateur… J’avais toujours admis que l’on avait à notre disposition un générateur de nombres aléatoires sur [0,1], qui génère à chaque appel des nombres indépendants, supposés uniformément distribués sur [0,1].

Je voulais prendre un peu de temps pour présenter quelques pistes de méthodes simples de méthodes permettant de générer des suites de nombre donnant l’illusion du hasard….
Une méthode plus ancienne (que j’ai connue, et oui) est celle des tables de nombres aléatoires, que l’on retrouvait par exemple dans les anciennes éditions du livre de Philippe Tassi ou de Gilbert Saporta. Je laisse les plus curieux traîner dans les vieux livres pour voir comment lire ce genre de table…. Une autre solution est celle que j’évoquais en introduction de mon (vieux) cours de séries temporelles à Dauphine (ici), à savoir l’idée de Slutsky d’utiliser les tirages du loto pour modéliser des prix (et donc inventer ainsi les processus MA, moving average), dans The Summation of Random Causes as the Source of Cyclical Processes. On peut trouver sur  le net les tirages du loto en France depuis des dizaines d’années (ici par exemple) et je suppose qu’en cherchant un peu, on peut trouver des choses similaires dans beaucoup d’autre pays.

Pour ma part, pour revenir sur une loi uniforme sur [0,1], je vais dans un premier temps diviser par 50.

Dans un second temps, on va concaténer deux boules consécutives (en base 50), puis se ramener sur [0,1].
Aussi, les boules 13 et 37 donnent comme nombre 0.2548.

On a ainsi un peu plus de possibilités sur les valeurs atteintes. Graphiquement, on obtient les résultats suivants,
XXXBref, on a à notre disposition des gens qui ont tiré des boules dans une urne. Alors rigoureusement, il faut faire attention car les tirages sont fait sans remise pour les 7 premières boules, puis on remet les boules pour le tirage suivant… Bref, les tirages ne sont pas rigoureusement indépendant les uns des autres (on ne peut pas avoir 3 fois consécutives le nombre “13”, alors qu’en théorie cela devrait être possible).
Bref, on est un peu limité avec ce genre d’outils, car ce qui m’intéresse c’est que mon ordinateur me génère autant de nombres aléatoires que je le souhaite.

Mexicoooooooooo

Après quelques jours de vacances (où contrairement à d’autres je n’ai pas emmené le Hamilton – cf ici – mais quelques très bons polars, dont le dernier Lehane, le dernier Indridason, ainsi que quelques vieux Rivages Noirs que j’avais ratés), je reviens en France, bien reposé…

Bon, j’en ai profité pour avancer deux papiers, reste à faire la mise en forme pendant les jours à venir…

The return period concept

In a paper with David Sibaï (that can be found here), we had a discussion about the concept of “return period” in hydrology. Actually, it looks like this concept has been introduced by Emil Gumbel in his book on Statistics of Extremes.

Graphs are also proposed,

The link between a probability and time can be established clearly using the geometric distribution. The time of the first success has the following distribution

\Pr(X = k) = (1 - p)^{k-1}\,p\,

where 0< p \leq 1 denotes the success probability. Then 
\mathrm{E}(X) = \frac{1}{p},
\qquad\mathrm{var}(X) = \frac{1-p}{p^2}.

This means that the time we have to wait – on average – before the first success is then simply the inverse of the success probability.
Note that this distribution is simply the discrete version of the exponential distribution, satisfying the memoryless property (in the context of continuous variates). And actually, this geometric distribution is the only discrete distribution satisfying such a property.
For instance, to illustrate this idea, the 1910 flood in Paris was suppose to be a centenial event. It does not mean that we must have a similar event next year, it means that we have to wait (still) 100 year – on average – before having a similar event. Or similarly, such an event occurs every year with a 1% probability, assuming temporal independence  (this was actually the point we discussed in our paper).

Looking for extremely rare events ?

As mentioned in a recent post on http://www.mafeco.fr/ (here), David Viniar who used to work as a director for Goldman Sachs mentioned that « We are seeing things that were 25-standard deviation events, several days in a row », (seehere). It was in 2007 much before the recent financial crisis (and the recent – real – extremal events)
I have to confess that I couldn’t do the maths with R, but this has been done recently in a research paper,
Dowd, K.; Cotter, J.; Humphrey, C.; Woods, M.(2008), « How unlucky is 25 sigma? », Journal of Portfolio Management, Vol.34(4), pp.76-80
It turns out that 25 times the standard deviation, for a gaussian distribution, is as unlikely yo happened as winning 21 or 22 consecutive weeks at the BritishNational Lottery ! So who can still believe that Gaussian distributions can still be used in financial econometrics ?

"sendo l'intento mio scrivere cosa utile a chi la intende…"