Tag Archives: computer

Removing Uncited References in a Tex File (with R)

Last week, with @3wen, we were working a the revised version of our work on smoothing densities of spatial processes (with edge correction). Usually, once you have revised the paper, some references were added, others were droped. But you need to spend some time, to check that all references are actually mentioned in the paper. For instance, consider the following compiled tex file :

Continue reading Removing Uncited References in a Tex File (with R)

La formation française, élitiste ?

En discutant avec @tomroud la semaine passée, j’ai découvert l’article Notre élite est toujours excellente, c’est le reste de la classe qui ne suit pas sur le site du Monde. C’est effectivement une idée reçue : si la formation est globalement aussi faible en France, c’est parce que la France serait un pays élitiste, et la formation de ses élites, en revanche, est une réussite. C’est – en gros – ce que raconte cet article. On pourra relire à ce sujet l’échange

Comment expliquer les résultats faibles en mathématiques, malgré la tradition d’excellence de la France dans ce domaine ?

Notre élite est toujours excellente, mais c’est le reste de la classe qui ne suit pas. Nos mauvaises notes à PISA ne nous empêcheront pas d’avoir des médailles Fields, mais elles nous empêchent d’avoir 100 % de citoyens efficaces dans leur travail et capables de se débrouiller dans un monde très mathématiques.

(dans la discussion, avec également, @Lazarrean@rdoume avait signalé un autre article dans Le Monde, la France championne des inégalités scolaires). Pourtant, ce n’est pas tout à fait ce que raconte le graphique, mentionné dans Surprising Test Results For Some Of The World’s Richest Students

où le “Ninth decile” (quantile à 90%) n’est pas si élevé que ça, comparé aux autres pays. Maintenant, pour raconter l’histoire complète, la discussion avec @tomroud est née de la lecture de La crise ne rend pas les riches plus riches et les pauvres plus pauvres, sauf en France,

En effet, @tomroud me disait qu’il avait vu le même genre de conclusion dans la comparaison des études PISA : en France, les bons sont meilleurs avec le temps, même si la moyenne, tout comme les niveaux des déciles les plus faibles, baisse. Dans une vision élitiste de l’éducation, ça ne serait pas gênant.

N’ayant rien trouvé me permettant de mieux comprendre , j’ai voulu aller regarder un peu plus en détails ce que disent vraiment les enquêtes PISA. Car elles sont disponibles en ligne (20122009200620032000). Bon, pour les lire plus rapidement, je vais les récupérer sur d’autres sites, pour avoir d’autres formats, plus simples à lire, en R,

> load(url("https://github.com/pbiecek/PISA2000lite/blob/master/data/math2000.rda"))
> load(url("http://beta.icm.edu.pl/PISAcontest/data/student2012.rda"))

Commençons par regarder un peu la distribution des notes, en mathématiques (je vais me restreindre aux tests de maths, qui sont souvent utilises pour dire que la France est effectivement un pays libelliste : les maths ont longtemps servi de critère, me semble-t-il).

> py="France"
> xf=sb[(sb$CNT==py)&(sb$ST04Q01=="Female")==TRUE,"PV1MATH"]
> xm=sb[(sb$CNT==py)&(sb$ST04Q01=="Male")==TRUE,"PV1MATH"]

Sur le graphique ci-dessous, on a la distribution des notes des filles (en haut à gauche) et des garçons (en bas à gauche). A droite, on a la fonction de survie : à une note donnée, on a en ordonnées la proportions d’élèves qui dépassent cette note. Par exemple 18,72% des garçons dépassent 600, et 12,42% des filles

> mean(xf>600)
[1] 0.1242105
> mean(xm>600)
[1] 0.1872207

Enfin, en bas à gauche, j’ai rajouté la fonction de moyenne en excédant (ou pour être plus précis, la Tail-Value-at-Risk), qui est la note moyenne que l’on peut espérer avoir si dans un certain décile. Par exemple, si on s’intéresse aux élèves dans le décile supérieur (dépassant le quantile à 90%), la note moyenne des garçons est 672 et celle des filles 649,

> mean(xm[xm>quantile(xm,.9)])
[1] 672.7608
> mean(xf[xf>quantile(xf,.9)])
[1] 649.281

Cette mesure est assez classique en tant que downside risk measure et permet de focaliser justement sur les queues de distribution (ce qui nous intéresse ici car on veut étudier les élites).

Si on zoom dans la partie à droite (les meilleures notes), on obtient

Mais il est difficile de faire des comparaisons avec d’autres pays sur ces graphiques. Par exemple, pour l’Australie, on obtient

 

Pour comparer les pays entre eux, on peut regarder la différence entre la note moyenne, et le niveau du quantile à 99% (on va chercher la vraie élite, les 1% des meilleurs… car c’est manifestement ce qui été évoque dans l’article, et qui intéresse du monde sur les inégalités de richesse)

L’élite en France est pas incroyablement meilleure que le cœur de la distribution, on en conviendra !

Si on compare un pays avec la France, on peut regarder la différence entre le quantile de niveau u (e.g. 50%, ou 90% ou 99%) entre ce pays et la France. En niveau absolu, ou en pourcentage.

> qpays("France",c(.5,.99))
      Female     Male
50% 496.8141 507.5635
99% 691.0734 717.2885

Sur les graphique ci-dessous, on a le différence absolue à gauche et relative à droite, avec l’Allemagne en haut, et l’Australie en bas,

On peut aussi regarder la Thaïlande et le Mexique (on est toujours en 2012)

Regardons un peu la Thaïlande par exemple. Les élèves (scolarises) médians sont 50 points en dessous des élèves Français. Mais pour le top 1%, les notes sont comparables… Donc manifestement, la formation des élites semble fonctionner dans d’autres pays que la France. En Australie, les inégalités scolaires semble encore plus fortes : les élèves médians sont globalement moins bon en Australie qu’en France (disons 15 points en dessous), mais l’élite est au dessus (de presque 50 points quand on regarde le top 1%).

Cela dit, comme j’avais dit qu’il pouvait aussi être intéressant de regarder les fonctions de moyennes en excédant (plus que les quantiles), on peut comparer ces fonctions,

Ici, les graphiques sont présentés en variation relative : au Mexique, la note moyenne d’une fille dans le top 20% (au dessus du quantile à 80%) est 17.5% plus faible qu’en France. On retrouve ici qu’en Australie par exemple, l’élite a – en moyenne – de meilleures notes qu’en France. Et en Thaïlande, plus on regarde les meilleurs élèves, plus les résultats sont comparables avec ceux des élèves français.

Qu’en est-il maintenant de l’évolution dans le temps ? Si on compare dans le temps, entre 2000 et 2012.

> sb2012=student2012[,c("PV1MATH","CNT","ST04Q01")]
> sb2000=math2000[,c("PV1MATH","CNT","ST03Q01")]
> sb2000=sb2000[!is.na(sb2000$PV1MATH),]
> sb2000$cnt=tolower(sb2000$CNT)
> sb2012$cnt=tolower(sb2012$CNT)
> ext2000=function(py){
+ x=sb2000[sb2000$cnt==py,"PV1MATH"]
+ v=c(mean(x),quantile(x,c(.1,.9)))
+ return(v)
+ }
> ext2012=function(py){
+ x=sb2012[sb2012$cnt==py,"PV1MATH"]
+ v=c(mean(x),quantile(x,c(.1,.9)))
+ return(v)
+ }
> listepays2000=unique(sb2000$cnt)
> listepays2012=unique(sb2012$cnt)
> R2000=sapply(listepays2000,ext2000)
> R2012=sapply(listepays2012,ext2012)
> D2000=data.frame(pays=listepays2000,
+ m2000=R2000[1,],
+ q12000=R2000[2,],
+ q92000=R2000[3,])
> D2012=data.frame(pays=listepays2012,
+ m2012=R2012[1,],
+ q12012=R2012[2,],
+ q92012=R2012[3,])
> D=merge(D2000,D2012)

Si on fait un dessin proche de celui sur les revenus (avec les évolutions des moyennes en bâtons, des 10% les plus bas en rond, et des 10% les plus hauts en rond), on obtient

Si effectivement le niveau moyen des étudiants baisse en France (comme dans plusieurs autres pays), on notera que l’idée que l’élite française s’en sort bien est un cliche sans réel fondement.

Computational Actuarial Science, with R

The book Computational Actuarial Science, with R is officially out. In the introduction of the book, and on the website of CRC, it is mentioned that the datasets can be found “in an R package on CRAN“, which is unfortunately incorrect. Some datasets are too large, so the package can not be uploaded on CRAN. Hopefully, Christophe host the package on his website.

> install.packages("CASdatasets", repos = "http://dutangc.free.fr/pub/RRepos/")

or

> install.packages("CASdatasets", repos = "http://dutangc.free.fr/pub/RRepos/", 
type = "source")

Here are the files :

Insurance datasets

A collection of datasets, originally for the book ‘Computational Actuarial Science with R’ edited by Arthur Charpentier (CAS with R). Now, the package contains a large variety of actuarial datasets.

Version: 0.9-8
Published: 2014-05-21
Author: Christophe Dutang
Maintainer: Christophe Dutang <christophe.dutang at ensimag.fr>
License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
NeedsCompilation: no

Downloads:

Reference manual: CASdatasets.pdf
Package source: CASdatasets_0.9-8.tar.gz
Package installation: go to this page
Windows binaries: r-release: CASdatasets_0.9-8.zip
OS X Snow Leopard binaries: r-release: CASdatasets_0.9-8.tgz
OS X Mavericks binaries: r-release: CASdatasets_0.9-8.tgz
Old sources: CASdatasets archive

Social Media Mining and Bioinformatics (with R)

In June and July, I receive copies of two books,

For the first one, two recent interesting books deal with the same topic. Reza Zafarani, Mohammad Ali Abbasi and Huan Liu published last year Social Media Mining: An Introduction. Actually, the book can be downloaded from dmml.asu.edu. And – of course – there is Matthew A. Russell‘s Mining the Social Web. The main interest of this new book seems to be that it should be just perfect for R users !

Continue reading Social Media Mining and Bioinformatics (with R)

Les anniversaires de vos amis sur Facebook

J’ai découvert avec un peu de retard le joli billet Les anniversaires de vos amis sur Facebook, qui tentait de répondre à la question

Si je possède n amis, quelle est la probabilité qu’il y ait au moins un jour dans l’année où je n’ai pas d’anniversaire à souhaiter ?

Ce problème, on peut aussi l’analyser en posant plutôt la question suivante,

Il y a en tout 365 jours dans une année. Combien d’amis doit-on espérer avoir pour qu’il n’y ait pas un jour sans fêter un anniversaire d’un de ses amis ?

Le problème est différent. Mais la résolution est assez jolie, si on regarde attentivement.

Si je note  le nombre d’amis qu’il me faut pour couvrir tous les anniversaires, et que je commence à regarder, ami après ami le nombre de jours couverts. Supposons que j’ai  jours de l’année couverts, et notons  le nombre de nouveaux amis qu’il va me falloir avant de couvrir  jours de l’année. On peut voir que

Le point intéressant avec cette notation, c’est qu’on connaît la loi de . En fait,  on reconnaît une loi géométrique

Aussi, on peut assez facilement avoir le nombre moyen d’amis qu’il va falloir pour couvrir tous les jours de l’année,

Je laisse les amateurs de calculs essayer de simplifier, mais ça me semble compliqué. Par contre, numériquement on obtient facilement

> sum(365/(1:365))
[1] 2364.646

En fait, la somme, on peut quand même la simplifier un peu, car on reconnaît

Et on sait que

où  est la constante d’Euler. Et effectivement, avec cette approximation, on obtient

> (log(365)+ 0.57721)*365
[1] 2364.144

Ce nombre d’amis est à rapprocher du calcul fait sur https://blogdemaths.wordpress.com/, qui obtenait

on a une chance sur deux pour qu’il y ait au moins un jour sans anniversaire à souhaiter, il faut avoir environ 2285 amis

Calculer un quantile est un peu plus compliqué pour moi (je ne peux calculer que des moyennes), mais on voit que les ordres de grandeurs sont assez proches.

Income distribution and Tour de France

A few days ago, Jean-François Mignot published an interesting article entitled Tour de France 2014 : pourquoi le vainqueur gagne 100 fois plus que le 10e. In this article, we have the following graph, with the income of the cyclist, as a function of his final ranking (the data where downloaded from http://sportbuzzbusiness.fr/)

> bike=read.csv(
+ "http://freakonometrics.free.fr/tourdefrance.csv",
+ sep=";",header=TRUE,dec=" ")

> bike[1:19,"prime"]=bike[1:19,"prime"]*1000
> plot(bike,log="y",type="b",
+ xlab="(Final) rank",ylab="Bonus")

As pointed out by Jean-François, if the winner gets a lot of money, the bonus decreases fast, very fast actually. Gini index is very high here

> library(ineq)
> ineq(X,type="Gini")
[1] 0.910461

and if we look at Lorenz curve, indeed, the Tour de France is not very equalitarian,

Continue reading Income distribution and Tour de France

London, Bayes and the Lloyd’s

Monday, we really had a great conference in London.

It was a great pleasure since I did learn a lot of things. And also a great honor to be the last speaker. Tuesday morning, I wanted to go to Thomas Bayes’grave, which is the the graveyard next to the CASS Business School. I had a good a apriori about where the grave should be,

but to be honest, it was not possible to get close enough to be able to read the name on it (even if I now know that it is the large one in the right lower corner of the picture)

Actually, on the internet, you can find some picture where the stone is clean, so you can learn that the grave is the “cotton” one – at least, you can easily read that name.

It was actualy more simple to see William Blake’s grave, as well as Daniel Defoe’s.

Then, with Leo, we’ve been to the Lloyd’s to see some friends, as well as Richard Rogers’s building.

At the 11th floor, you have a lot of rooms for meetings, as well as old paintings, to tell a bit more about the history of the company,

The building is just amazing. Unfortunately, to get in, there is a dress code. A sort of strict one actually. Leo is working for RBC, so he casually wears a suit. But I don’t. I mean, I did have a shirt, but as someone mentioned, “there is no collar !” (I don’t want to put my friend into trouble for helping me getting in).

So, after going throught the basement, we’ve been able to reach the elevator, and go on top.

The building is not exactly located where Edward Lloyd got his coffee shop (even after moving at the end of 1691 on Lombard street), but the Lloyd’s is still a legend for anyone interested in the history of insurance, and more generally, the history of risk modeling (and management).

Computational Actuarial Science

After some delay, the book Computational Actuarial Science with R is now annonced for July 2014. I don’t know if we will be able to get copies for the R in Insurance conference, in London, but I guess everyone is working on it. And kindly, CRC sent me the following flyer, with some reduction code to save 20% (when ordering on CRC’s website).

Modèles de Prévision

Cet session d’été, je redonne le cours ACT6420, modèles de prévision. Le plan de cours est maintenant en ligne (comme souvent, les salles de cours ne sont pas encore attribuées définitivement mais elles seront disponibles sur le site de l’UQAM). Les transparents pour la première partie de cours sont en ligne (il s’agit de rappels de statistique, et d’une introduction aux modèles de régression).

Des exercices et des codes informatiques seront mis en ligne sur /courses/modeles-de-prevision.

Durant le cours, j’essayerais d’illustrer avec le maximum d’exemples sur des données réelles, à partir de codes R. Quelques séances de démonstrations se feront en salle informatique, afin de découvrir les rudiments de R. Pour aller plus loin, je recommande

  • “R pour les débutants”d’Emmanuel Paradis (PDF)
  • “Brise Glace-R” d’Andrew Robinson er Arnaud Schloesing (PDF)
  • “Introduction à la programmation en R” de Vincent Goulet (PDF)

ainsi que quelques tutoriels en ligne, par exemple, dans un ancien billet, mais aussi sur https://youtube.com/.

Copula Density Estimation

The joint paper, written with Gery Geenens and Davy Paindaveine, entitled Probit transformation for nonparametric kernel estimation of the copula density” is now online on http://arxiv.org/abs/1404.4414

Copula modelling has become ubiquitous in modern statistics. Here, the problem of nonparametrically estimating a copula density is addressed. Arguably the most popular nonparametric density estimator, the kernel estimator is not suitable for the unit-square-supported copula densities, mainly because it is heavily affected by boundary bias issues. In addition, most common copulas admit unbounded densities, and kernel methods are not consistent in that case. In this paper, a kernel-type copula density estimator is proposed. It is based on the idea of transforming the uniform marginals of the copula density into normal distributions via the probit function, estimating the density in the transformed domain, which can be accomplished without boundary problems, and obtaining an estimate of the copula density through back-transformation. Although natural, a raw application of this procedure was, however, seen not to perform very well in the earlier literature. Here, it is shown that, if combined with local likelihood density estimation methods, the idea yields very good and easy to implement estimators, fixing boundary issues in a natural way and able to cope with unbounded copula densities. The asymptotic properties of the suggested estimators are derived, and a practical way of selecting the crucially important smoothing parameters is devised. Finally, extensive simulation studies and a real data analysis evidence their excellent performance compared to their main competitors.”

Computational Actuarial Science

Last week, we’ve been through the book, completely, one last time, before sending it back to the publisher, with some comments and remarks, before publication ! So, this is it, the book will finally appear soon ! It was schedule for this week actually, but… you know. It should appear sometime by the end of May, or beginning of June. I will keep you posted on this blog.

A few months ago, we published with Christophe Dutang an ebook on the same topic, in French, online on cran.r-project.org/doc/contrib/. This contribution was based on lecture notes we had. When I’ve been asked to publish an English version, by John Kimmel, I was honored, but I thought it would be some kind of fraud if I write a book on that topic. I do know a bit of actuarial science, and a bit of R, but most of the advanced computations rely on others packages. Because I am extremely lazy, I did not try (so far) to edit my own package. I frequently publish some lines of codes on my blog, but nothing too serious.

So, for this book, I decided to ask those who actually did publish a package used in actuarial computations (or who did work previously on packages comparison for instance) to write a chapter, in this book. I am usually not a big fan of books with twenty contributors, because there is no coherence. So here, my task was to link all those chapters together, to make sure that notations are coherent, etc. Over 700 pages, that was difficult. And I asked all of them not only to illustrate actuarial concepts with some R code, but also give – if possible – some self written function, to understand the algorithm, but also some built-in functions. The goal was to explain the core of the algorithm. Some codes might not be efficient, but they help to understand how it could be possible to compute some actuarial quantities.

The first chapter is an Introduction (to the R language) I wrote with Rob Kaas (everyone in the actuarial community knows Rob ! not only as the Editor of Insurance: Mathematics & Economics. but also as a prolific author, including the popular textbook Modern Actuarial Risk Theory – Using R). The aim is to help those who might use another language for actuarial computation to understand the basics of the grammar, to read and write in R. I will probably publish a longer post, to explain the structure of that chapter. And show some codes.

  • Methodology 

The first section is a very general methodology section. It starts with Standard Statistical Inference by Christophe Dutang (Christophe is extremely active in the R community, as the maintainer of the Distributions task view page, for instance).Then, Ben Escoto (Ben works in the insurance industry, and launched the actuarial vignettes in R a few years ago) and myself, wrote a chapter which can be seen as an introduction to the Bayesian Philosophy for actuaries (I will give a talk on that topic at the R in Insurance conference this summer, so additional material will be online soon). With Stéphane Tufféry, we Statistical Learning (Stéphane published a Data Mining and Statistics for Decision Making a few years ago). Then, I wanted a chapter dedicated to Spatial Analysis. I did ask Renato Assunção, Marcelo Azevedo Costa, Marcos Oliveira Prates, and Luís Gustavo Silva e Silva to write that chapter (I met them a few years ago while I was visiting Renato in Belo Horizonte, while that started to work on spatial aspects of actuarial science). And finally, Eric Gilleland and Mathieu Ribatet wrote the chapter on Reinsurance and Extremal Events (both of them work on climate and extreme values, and they did publish a very interesting software review for extreme value analysis a few years ago).

  • Life Insurance

For the section on life insurance, I asked Giorgio Spedicato to write the chapter on Life Contingencies (Giorgio is the author of the lifecontingencies package). Then Heather Booth, Rob J. Hyndman, and Leonie Tickle agreed to write the chapter on Prospective Life Tables (here we have a great match, with a demographer, an actuary, and… Rob. Every one who studied time series knows Rob. He is the author of the amazing forecast package, as well as the demography package, among many others. And he has a great blog too). To go further, Julien Tomas and Frédéric Planchet wrote the chapter on Prospective Mortality Tables and Portfolio Experience (both of them published the ELT – experience life tables – package a few months ago). And finally, there is a chapter on Survival Analysis by Frédéric Planchet  and Pierre-E. Thérond (they did publish a book – in French – on survival analysis for actuarial science, with examples in R).

  • Finance

For the section on financial computations, Yohan Chalabi and Diethelm Würtz wrote two chapters, one on Stock Prices and Time Series and one on Portfolio Allocation (both of them have worked on the Rmetrics project, with the timeSeriesfArmafGarchfPortfolio packages). And Sergio S. Guirreri wrote a chapter on Yield Curves and Interest Rates Models (Sergio is the author of the YieldCurve package).

  • Non-Life Insurance

Last, but not least, there is a section on non-life insuranceJean-Philippe Boucher (who published several articles on counts models) and myself, wrote the chapter on General Insurance Pricing. Then, I asked Katrien Antonio, Peng Shi, and Frank van Berkum to go further, with a chapter on Longitudinal Data and Experience Rating (I know Katrien from my PhD, and she was already working on that topic by that time… she did publish great surveys on that topic). And finally, Claims Reserving and IBNR is a chapter I wanted to write, because it’s a topic I love, but I asked Markus Gesmann to write it (Markus is known not only for his googleVis package, but also for the ChainLadder package – not to mention his awesome blog).

I will try to post some additional material on this blog, with R code (of course), graphs, and slides. And probably some pdfs with answers for the exercises. And all the datasets will be available in a CASdatasets package (online soon).

On retourne au pôle nord ?

Depuis plusieurs jours, j’essaye de me préparer autant que possible au débat de mardi soir, et plus j’avance dans mes lectures, plus je trouve que la question “à qui appartient le pôle nord” est complexe. Un autre exemple récent de “revendication” d’une région par d’autres pourrait être celui de l’Ukraine. Pour simplifier outrageusement un débat (là aussi), j’ai l’impression qu’une partie de l’Ukraine se sent européenne, et une autre se sent russe. Un outil pour visualiser cette appartenance peut être d’utiliser des données de langue, par exemple.

File:Ukraine cencus 2001 Ukrainian.svg

Avec des données et une visualisation sur une carte, on peut tenter d’éclairer un peu. Mais dans le cas du pôle nord, on ne peut pas vraiment utiliser de données de ce genre (pas à ma connaissance, malheureusement). Cela dit, on peut apprendre des choses passionnantes dans un ancien numéro de Courrier International sur cette même question :

Dans mon précédant billet, j’avais parlé de la discussion que j’avais eu, rapidement, avec ma fille. En emmenant mon fils à l’escrime vendredi soir, je lui expliquais la problématique du débat, et sa réponse a été sans ambiguïté: le pôle nord appartient à celui qui a mis son premier un drapeau dessus ! On retrouve là la visualisation graphique de l’affiche du débat. Maintenant je dois avouer que ce genre de revendication me gêne un peu. La symbolique de la pose du drapeau est troublante quand on y pense. Une des images à laquelle on pense immédiatement est d’ailleurs le drapeau américain, posé sur la lune (en écho au drapeau hissé lors de la bataille d’Iwo Jima, 硫黄島, en février et mars 1945). A l’époque du drapeau américain sur la lune, il avait d’ailleurs été précisé « this act is intended as a symbolic gesture of national pride in achievement and is not to be construed as a declaration of national appropriation by claim of sovereignty » dans un texte adopté par le Sénat. Donc peut être que cette réponse est trop rapide ! Je voulais donc revenir sur deux éléments clés du débat : les pays qui pourraient revendiquer le pôle nord, et le pôle nord, en tant que tel.

  • les pays autour du pôle nord

Pour l’instant, j’avais fait quelques cartes pour visualiser le pôle, par rapport aux pays qui l’entourent (essentiellement le Canada, la Russie et le Groenland). J’avais retenu comme contour pour mes cartes les limites terrestres des pays, tels que définis sur les cartes. Mais on peut aller plus loin, comme l’évoquent certains, en tenant compte de la zone de 200 milles maritimes (la notion d’Exclusive Economic Zone).

La (petite) région au centre est celle qui est découpée dans la carte de Courrier International, en début de billet. Comme on le voit sur la carte ci-dessus, en rajoutant ces zones, effectivement, nos trois nations sont très proches du pôle (le contour de glace est celui observé pendant les 12 mois de 2013),

Pour moi, l’avantage de ces zones est que c’est plus simple de regarder l’intersection entre la zone de glace, et ces zones économiques (et de calculer l’aire) que de mesurer la longueur de frontière entre un pays et la zone de glace. On va essayer en fin de billet de regarder quelle proportion de la zone de glace chevauchait une zone exclusive (je passe ici sous silence la difficulté technique de calcul d’une aire d’une région contenant le pôle nord… étrangement, c’est complexe, d’un point de vue computationel. Les plus curieux iront voir le billet dédié).

Cela dit, certains pourraient critiquer l’idée de tenir compte de cette extension (horizontale) du territoire, en tenant compte de ces zones économiques. C’est quoi la prochaine étape ? Une appropriation verticale, afin de prendre possession des nuages. Cela dit, je plaisante à peine car la possession verticale de l’espace existe, avec les espaces aériens… D’ailleurs les russes ont évoqués une appropriation verticale, mais sous terre, car la dorsale de Lomonossov  – Хребет Ломоносова – relirait le pôle directement au territoire russe.

Mais je vais arrêter là ma discussion sur les pays.

  • quel pôle nord ?

J’avais un peu laissé de côté la question, mais je pense qu’il serait quand même intéressant d’en parler au moins 2 minutes. De quoi parle-t-on quand on parle du pôle nord ? Pour l’instant, j’avais utilisé la définition usuelle du pôle nord, correspondant au point le plus septentrional de la terre (ce qui, j’en conviens, est un peu tautologique comme définition). En terme de coordonnées, il s’agit du point à 90° de latitude Nord. Il n’a pas de longitude, car c’est précisément le point au nord où tous les méridiens (et tous les fuseaux horaires) se coupent. J’avais aussi dit, dans mon précédant billet qu’il s’agissait du point de la surface terrestre, de l’hémisphère nord, situé sur l’axe de rotation de la terre. Ce point géographique n’est malheureusement pas fixe à la surface de la Terre, car l’axe de rotation de la terre varie (très faiblement).

Mais il existe un autre pôle, le pôle nord magnétique, qui est l’unique sur la surface de la terre où le champ magnétique (terrestre) pointe vers le centre de la terre (et orthogonal à la surface de la terre, s’il pointe vers le haut, c’est le pôle sud). A priori, on pourrait se dire que ce point n’aura pas grand intérêt dans la discussion, car s’approprier un point mouvant serait surprenant. Car il bouge beaucoup, en témoigne la carte ci-dessous (les coordonnées sont facilement accessibles, sur http://ngdc.noaa.gov/geomag/data/poles/ par exemple)

Un utilisant ces champs, on peut d’ailleurs définir un pôle nord gravitationnel (je me permets d’inventer le mot, par analogie au pôle nord magnétique), à l’aide non pas du champ magnétique, mais du champ gravitationnel. En anglais, on parle de polar wander. Pour faire simple, ce pôle nord sera lié à la répartition de la masse, sur terre, et des mouvements d’eau, en particulier (le reste étant assez stable à côté). Or comme l’ont montré récemment Chen, Wilson, Ries et Tapley, dans Rapid ice melting drives Earth’s pole to the east (article remarqué par Richard Lovett, en mai 2013, dans Nature) le pôle nord semble se déplace plus rapidement que prévu, à cause du réchauffement climatique.

Mais je pense que cette discussion (quoi que potentiellement passionnante) nous éloigne de ce qui devrait être évoqué mardi soir. En fait, mon point serait que, plutôt que de visualiser le pôle nord, en tant que point, il est plus intéressant de visualiser une région. Et plutôt que de considérer une région au delà du cercle polaire (j’avais utilisé des cercles dans mon précédant billet), on peut davantage utiliser la zone glacée, comme le montre la carte

Encore une fois, le soucis est que cette zone glacée varie énormément ! On observe qu’entre la moitié et les deux-tiers de la surface de la glace se situe dans les zones économiques des trois principaux pays limitrophes, la Russie et le Canada.

Et si on y regarde de plus près, on observe d’ailleurs que les surfaces occupées respectivement par la Russie et le Canada sont comparables,

avec toutefois un comportement cyclique, la présence canadienne étant au plus faible à l’automne,

J’ai l’impression que je pourrais passer des heures sur ces cartes. Mais je vais probablement m’arrêter là pour ma préparation du débat. Et même si ces débats ne laissent que peu de temps pour la prise de notes, j’essayerais de prendre un peu de temps pour revenir sur ce qui s’est dit, en milieu de semaine, sur le blog, ou sur Twitter.

Temperatures Series as Random Walks

Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl)
library(XML)
options(RCurlOptions = list(useragent = "R"))
HEURE=0:23
extracttemp=function(Y,M,D){
url=paste(
"http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=",
M,"&Day=",D,sep="")
wp <- getURLContent(url)
doc <- htmlParse(wp, asText = TRUE) 
docName(doc) <- url
tmp <- readHTMLTable(doc)
basejour=data.frame(Year=Y,Month=M,Day=D,
Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25])
return(basejour)}
B=NULL
for(y in 1955:2013){
for(d in 1:31){
B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)")
I=which(B$Year==2013)
lines(B$X[I],B$Temp[I],col="red")

In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing.

I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done)

correcty=function(Y){
I=which(is.na(Y))
	if(length(I)==0){Yc=Y}
	if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]}	
return(Yc)
}

Now, we can compute the p-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary)

DF=matrix(NA,2013-1954,3)
library(urca)
for(y in 1955:2013){
Z=B$Temp[which(B$Year==y)]
	Zc=correcty(Z)
	DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value)	
	DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value)	
	DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value)
}

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF
DFP[,1]=DF[,1]<.05
DFP[,2:3]=DF[,2:3]>.05
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")
plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="")
axis(1)
text(1952,.5,"KPSS")
text(1952,1.5,"PP")
text(1952,2.5,"ADF")
for(y in 1955:2013){
for(i in 1:3){
polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).