# 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.

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.

# R in Insurance, July 2014, London

As mentioned a few months ago, the second conference on R in Insurance will be held on Monday 14 July 2014 at Cass Business School in London, UK. Registration is open.

# 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
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")
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).

# Le Vote par Procuration en France

La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.

Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.

• Modèle de régression, et recherche de variables explications

Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que $P_i\sim\mathcal{B}(N_i,p_i)$ où $P_i$ est le nombre de procurations dans le bureau de vote $i$, et où $N_i$ est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que $p_i$, la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT > baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

Si on suppose que $p_i$ est fonction $X_i$ le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference] > baset2$taux=baset2[,variable]/baset2[,reference]

il est légitime de tenter une régression logistique,

$p_i=h(X_i)=\frac{\exp[\beta_0+\beta_1 X_i]}{1+\exp[\beta_0+\beta_1 X_i]}$

voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,

$p_i=\tilde h(X_i)=\frac{\exp[s(X_i)]}{1+\exp[s(X_i)]}$

Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)

> b=hist(baset1$taux,plot=FALSE) > library(splines) > regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1) > regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2) > u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001) > ND=data.frame(taux=u) > ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){ + polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)} > polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)), + border=NA,density=30,col=CL[4]) et, pour la régression logistique standard (linéaire) > lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2) > regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1) > regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2) > ND=data.frame(taux=ug) > pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response") > pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response") > lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2) > legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)], + lwd=2,lty=1,border=NA) (en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration). On peut faire la même chose sur le taux de logement HLM, dans le quartier, Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de $\beta_1$ laisse à désirer (toute suggestion est la bienvenue !). > summary(regt1l) Call: glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, data = baset1, weights = INSCRITS) Deviance Residuals: Min 1Q Median 3Q Max -12.9549 -1.5722 0.0319 1.6292 13.1303 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.70811 0.01516 -244.6 <2e-16 *** taux 1.49666 0.04012 37.3 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 12507 on 836 degrees of freedom Residual deviance: 11065 on 835 degrees of freedom AIC: 15699 Number of Fisher Scoring iterations: 4 > summary(regt2l) Call: glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, data = baset2, weights = INSCRITS) Deviance Residuals: Min 1Q Median 3Q Max -15.4872 -1.7817 -0.1615 1.6035 12.5596 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.24272 0.01230 -263.61 <2e-16 *** taux 1.45816 0.03266 44.65 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 9424.7 on 836 degrees of freedom Residual deviance: 7362.3 on 835 degrees of freedom AIC: 12531 Number of Fisher Scoring iterations: 4 On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire $\frac{P_i}{N_i}=\beta_0+\beta_1 X_i+\varepsilon_i$(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici > regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1) > regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2) > u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001) > ND=data.frame(taux=u) > ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){ + polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)} > polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)), + border=NA,density=30,col=CL[4]) > lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2) > regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1) > regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2) > ND=data.frame(taux=ug) > pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response") > pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response") > lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2) > legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)], + lwd=2,lty=1,border=NA) (j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre) Cette fois, on a une sortie de régression plus classique, > summary(regt1l) Call: lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS) Weighted Residuals: Min 1Q Median 3Q Max -1.9994 -0.2926 0.0011 0.3173 3.2072 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.021268 0.001739 12.23 <2e-16 *** taux 0.054371 0.004812 11.30 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.646 on 835 degrees of freedom Multiple R-squared: 0.1326, Adjusted R-squared: 0.1316 F-statistic: 127.7 on 1 and 835 DF, p-value: < 2.2e-16 > summary(regt2l) Call: lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS) Weighted Residuals: Min 1Q Median 3Q Max -2.9029 -0.4148 -0.0338 0.4029 3.4907 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.033909 0.001866 18.17 <2e-16 *** taux 0.079749 0.005165 15.44 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.6934 on 835 degrees of freedom Multiple R-squared: 0.2221, Adjusted R-squared: 0.2212 F-statistic: 238.4 on 1 and 835 DF, p-value: < 2.2e-16 On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions). • Visualiser les bureaux de vote à Paris Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net) > library(maptools) > library(rgdal) > library(classInt) > paris=readShapeSpatial("paris-cartelec.shp") Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes > elec=data.frame() > elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1]
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

Reste à conclure, en faisant une visualisation graphique de nos données

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")),
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

Dans ce cas, la visualisation des résidus donne

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1]
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")),
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.

# Personal Analytics with RSS Feeds

I am currently working on a paper on Academic Blogging, from my own experience. And I wanted to do something similar to Stephen Wolfram’s personal analytics of my life. More specifically, I wanted to understand when I do post my blog entries. If I post more entries during office hours, then it should mean that, indeed, I consider my blog as a part of my job (which is something I believe, actually). On the other hand, if I post more in the evening, or in the middle the night, then it could mean that my blog is clearly only for fun, and somehow outside the official academic time schedule.

With the help of @3wen, we have here a function that can read rss feeds, and extract the publication date (and other pieces of information actually),

> library(XML)
> library(dplyr)
+   if(length(doc)>1){
+   lesArticles <- xpathApply(r <- xmlRoot(doc), "//item")
+   infosUneEntree <- function(x){
+   title <- sapply(xpathApply(x, "//title"), xmlValue)
+   pubDate <- sapply(xpathApply(x, "//pubDate"), xmlValue)
+ }
+ df <- lapply(lesArticles, infosUneEntree)
+ df <- data.frame(do.call("rbind", df))
+ return(df)
+ }
+ else{return(NA)}
+ }

The trick is that the page containing the rss feeds is truncated: you get only 30 post (the latest ones). With WordPress, you can easily go further (thanks @3wen) using

> df.freak2 <- baseRSS("http://freakonometrics.hypotheses.org/feed?paged=2")
Namespace prefix dc on creator is not defined
Namespace prefix content on encoded is not defined
Namespace prefix wfw on commentRss is not defined
Namespace prefix slash on comments is not defined
> head(df.freak2)
title
1       S\303\251ries chronologiques, syllabus
2 Copules et valeurs extr\303\252mes, syllabus
3         Jimmy, Mile End, et le Qu\303\251bec
4                Multivariate Archimax copulas
5                     Somewhere else, part 107
6     Informatique (sans ordinateur), partie 1
1 http://freakonometrics.hypotheses.org/11593 Mon, 06 Jan 2014 00:31:52 +0000
2 http://freakonometrics.hypotheses.org/11595 Mon, 06 Jan 2014 00:31:21 +0000
3 http://freakonometrics.hypotheses.org/11362 Sun, 05 Jan 2014 03:33:31 +0000
4  http://freakonometrics.hypotheses.org/7673 Sat, 04 Jan 2014 11:01:05 +0000
5 http://freakonometrics.hypotheses.org/11584 Fri, 03 Jan 2014 15:34:29 +0000
6 http://freakonometrics.hypotheses.org/11138 Fri, 03 Jan 2014 07:15:03 +0000

and if we try to get a page that does not exist, we got the following error

> df.freakFaux <- baseRSS("http://freakonometrics.hypotheses.org/feed?paged=2000")
Error : 1: failed to load HTTP resource

(unfortunately, I could not do it with https://feeds.feedburner.com/ for instance). With the following code, we can extract information about all the posts online on my blog

> df.freak <- NULL
> for(i in 1:2000){
+   df.tmp <- baseRSS(paste("http://freakonometrics.hypotheses.org/feed?paged=", i, sep = ""))
+   if(length(df.tmp)>1){
+     df.freak <- rbind(df.freak, df.tmp)
+   }else{ break }
+ }

All that is just fine. Now, let us write a small function to convert the date into some format I can use (here, I want to study the hour, as well as the week day).

> LD=c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")
> datahour=function(txt){
+ wd=substr(as.character(txt),1,3)
+ wdy=which(LD==wd)
+ y=substr(as.character(txt),13,16)
+ h=substr(as.character(txt),18,19)
+ mn=substr(as.character(txt),21,22)
+ T=as.numeric(h)+as.numeric(mn)/60
+ return(data.frame(weekday=wdy,time=T,year=as.numeric(y)))}

+ L=unlist(lapply(as.character(df$pubDate),datahour)) + db=data.frame( + D=L[names(L)=="weekday"], + T=L[names(L)=="time"], + Y=L[names(L)=="year"]) + return(db)} Here, I extract the week day, the time in the day (continuous, from 0 till 24, excluded). With the following function we can see the proportion of posts per week day, > hc=rev(heat.colors(100)) > weekday=function(db,yearinf=FALSE){ + y=unique(db$Y)
+ if(yearinf==TRUE) y=y[-which.max(y)]
+ if(yearinf==FALSE) y=y[-c(which.max(y),which.min(y))]
+ L=NULL
+ for(i in y){
+ sB=subset(db,db$Y==i) + L=rbind(L,table(sB$D)/nrow(sB)*100)}
+ barplot(t(L[nrow(L):1,]),names=rev(y),col=hc[c(rep(15,5),rep(70,2))])
+ }

(from the bottom to the top, Monday till Friday in light yellow, and Saturday and Sunday in light red). Here, on my own blog, it would be

> weekday(datarss(df.freak))

For the hour, it was slightly more technical (I could not find a decent and simple way to plot the graph I was looking for graph, so I did it by myself)

> hour=function(db,yearinf=FALSE){
+ y=unique(db$Y) + if(yearinf==TRUE) y=y[-which.max(y)] + if(yearinf==FALSE) y=y[-c(which.max(y),which.min(y))] + L=NULL + for(i in y){ + sB=subset(db,db$Y==i)
+ if(i==2013) t=table(floor((sB$T+6)%%24))/nrow(sB)*100 + if(i<2013) t=table(floor(sB$T))/nrow(sB)*100
+ t=t[as.character(0:23)]
+ names(t)=as.character(0:23)
+ t[is.na(t)]=0
+ L=rbind(L,t)}
+ plot(y,rep(24,length(y)),ylim=c(-3,24),axes=FALSE,
+ xlim=c(min(y)-.5,max(y)+.5),xlab="",ylab="",col="white")
+ axis(2)
+ for(i in y){
+ text(i,-2,i)
+ for(j in 0:23){
+ polygon(c(i-.4,i-.4,i+.4,i+.4),
+ c(j,j+1,j+1,j),border=NA,col=hc[L[max(y)-i+1,j+1]/max(L)*98+1])
+ }}}

Just a short comment here. If you look at the code, there is a difference between 2013, and before. The reason is simple: in December 2012, I officially decided to migrate from my old blog to this new one. All the post prior December 2012 were initially published on the old blog. Which was at Montréal (East Coast) time. And I have the feeling that my new blog has a European time. So I did translate, of 6 hours. But the problem might be more complicated actually

> hour(datarss(df.freak))

Now, if we try to comment. On the week days, I find it a bit scary, to see that I spend so much time during the weekends on my (supposed to be) professional blog. And on the hour, I can explain the 2013 easily. I usually spend most of my evenings working (on the blog, or on my courses, or on my research). But usually, I try to avoid posting an entry at 2 a.m. So usually, I keep it until the morning, then when I arrive at the office, I finalize the post, and I make it available.

To understand the difference with previous years, I should probably add a technical comment : the previous blog was on a dotclear platform. On dotclear the Publication time is not exactly the time the post was officially posted online, but the default value is more the time the post was saved for the first time. So there might be some slight differences. I believe that previously, I started to work on a post in the afternoon, then I might spend some time in the evening, even the day after, but when I publish it, if I do not change the default settings, then the publication time would be the afternoon, when I did save the post.

Let us try on another blog… The problem is that is it is quite difficult to get old entries from the rss feeds. Except with WordPress… So I tried to run the previous code on http://economix.blogs.nytimes.com/. The extraction is simple here.

 501 Tue, 14 May 2013 04:01:

But here again, I do have trouble with 2013. To be more specific, when I look at the feeds I get

while I have on my side

 501 Tue, 14 May 2013 04:01:29 +0000
502 Mon, 13 May 2013 19:51:06 +0000
503 Mon, 13 May 2013 04:01:46 +0000
504 Fri, 10 May 2013 21:11:58 +0000
505 Fri, 10 May 2013 18:45:48 +0000
506 Fri, 10 May 2013 17:33:55 +0000
507 Fri, 10 May 2013 13:00:41 +0000
508 Fri, 10 May 2013 04:01:53 +0000

I have here a 4 hours difference I cannot explain. But it looks fine before 2013. If I use the previous code, with (in the loop)

+   df.tmp <- baseRSS(paste("http://economix.blogs.nytimes.com/feed/?paged=", i, sep = ""))

we can get, for instance the following  graph,

We do observe an interesting dynamics here : I guess that previously people were working during the day, and then posting at the end of the day. It looks like, now, people work in the day, sometimes late in the evening, but wait till the next morning to post the entry. Just as I did, in order to read one last time, with a fesh mind… Anyway, I still have to understand what did happened in 2013, just to make sure that the data I extract can be used…

# Inference for MA(q) Time Series

Yesterday, we’ve seen how inference for $AR(p)$ time series was possible.  I started  with that one because it is actually the simple case. For instance, we can use ordinary least squares. There might be some possible bias (see e.g. White (1961)), but asymptotically, estimators are fine (consistent, with asymptotic normality). But when the noise is (auto)correlated, then it is more complex. So, consider here some $MA(2)$ time series

$X_t=\varepsilon_t +\theta_1 \varepsilon_{t-1}+\theta_2\varepsilon_{t-2}$

for some white noise $(\varepsilon_t)$.

> theta1=.25
> theta2=.7
> n=1000
> set.seed(1)
> e=rnorm(n)
> Z=rep(0,n)
> for(t in 3:n) Z[t]=e[t]+theta1*e[t-1]+theta2*e[t-2]
> Z=Z[800:1000]
> plot(Z,type="l")

• Using the empirical autocorrelations

The first idea might be to use the first two (empirical) autocorrelations (the two that are supposed to be – theoretically – non null).

$\rho(1)=\frac{\theta_1+\theta_1\theta_2}{1+\theta_1^2+\theta_2^2}$

$\rho(2)=\frac{\theta_2}{1+\theta_1^2+\theta_2^2}$

with $\rho(h)=0$ when $h\geq0$. We also have the following relationship on the variance of the process

$\text{Var}(X_t)=\sigma^2\cdot(1+\theta_1^2+\theta_2^2)$

With those three equations, for three unknown parameters, $\theta_1$$\theta_2$ and $\sigma$, we simply have to solve (numerically) that system of equations,

> v=c(as.numeric(acf(Z)$acf[2:3]),var(Z)) > v [1] 0.1658760 0.3823053 1.6379498 > library(rootSolve) > seteq=function(x){ + F1=v[1]-(x[1]+x[1]*x[2])/(1+x[1]^2+x[2]^2) + F2=v[2]-(x[2])/(1+x[1]^2+x[2]^2) + F3=v[3]-(1+x[1]^2+x[2]^2)*x[3]^2 + return(c(F1,F2,F3))} > multiroot(f=seteq,start=c(.1,.1,1))$root
[1] 0.1400579 0.4766699 1.1461636

$f.root [1] 7.876355e-10 4.188458e-09 -2.839977e-09$iter
[1] 5

$estim.precis [1] 2.605357e-09 We are a bit far away from the true values, used to generate our sample. And if we consider 1,000 sample (instead of only one), we still have the bias, and a large variance for our three estimators, • Using least square techniques We can try something quite different here. The problem we have is that we do not observe the noise $(\varepsilon_t)$, we only observe our series $(X_t)$. But we can try to rebuild that series (call it $(u_t)$ since we’re not sure it will be a reconstruction of the noise). As suggested in Box & Jenkins (1967), assume that the first two values are null. And then, use $u_t =X_t-\left(\theta_1u_{t-1}+\theta_2un_{t-2}\right)$ and then, we can use least square techniques $\min_{(\theta_1,\theta_2)}\left\{\sum_{i=1}^tu_i^2\right\}$ The code will be > V=function(p){ + theta1=p[1] + theta2=p[2] + u=rep(0,length(Z)) + for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2] + return(sum(u^2)) + } If we try to minimize the sum of the squares of the residuals, we get > optim(par=c(.1,.1),V)$par
[1] 0.2751667 0.6723909

$value [1] 225.8104$counts
77       NA

$convergence [1] 0$message
NULL

which is close to the true value. Another good thing is that, if we compare that rebuilt noise with the true one (since we actually have it), then we have the same vector,

> plot(e[800:1000],col="blue",type="l")
> theta1=0.2751667
> theta2=0.6723909
> u=rep(0,length(Z))
> for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2]
> lines(1:201,u,col="red")

So far, so good. And if we look at 1,000 samples, we get

It looks like we have some bias here. And since the two estimators should be negatively correlated, one over-estimates, while the other one under-estimates.

• Using the (global) maximum likelihood technique

And a final method might be to use the maximum likelihood technique (globally). Again, if we assume that we have a Gaussian i.i.d noise, then the vector $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ is Gaussian, with a simple variance matrix (since a lot of elements will be null),

> library(mnormt)
> GlobalLogLik=function(A,TS){
+ n=length(TS)
+ theta1=A[1];  theta2=A[2]
+ sigma=A[3]
+ SIG=matrix(0,n,n)
+ rho=rep(0,n)
+ rho[1]=1
+ rho[2]=(theta1+theta1*theta2)/(1+theta1^2+theta2^2)
+ rho[3]=(theta2)/(1+theta1^2+theta2^2)
+ for(i in 1:n){for(j in 1:n){
+ SIG[i,j]=rho[abs(i-j)+1]}}
+ gamma0=(1+theta1^2+theta2^2)*sigma^2
+ SIG=gamma0*SIG
+ return(dmnorm(TS,rep(0,n),SIG,log=TRUE))}
> LogL=function(A) -GlobalLogLik(A,TS=Z)
> optim(c(.1,.1,1),LogL)
$par [1] 0.2584144 0.6826530 1.0669820$value
[1] 298.8699

$counts function gradient 86 NA$convergence
[1] 0

$message NULL Here, the values that minimize the likelihood are rather close to the ones used to generate our sample. And if we run this algorithm on 1,000 samples, we can see that those estimates are fine, I could not find other ideas, to estimate those parameters. I guess we can use the partial autocorrelation function, since we have relationships that can be related to Yule-Walker equations for $AR(p)$ time series. # Inference for AR(p) Time Series Consider a (stationary) autoregressive process, say of order 2, $Y_t =\varphi_1 Y_{t-1}+\varphi_2 Y_{t-2}+\varepsilon_t$ for some white noise with variance . Here is a code to generate such a process, > phi1=.25 > phi2=.7 > n=1000 > set.seed(1) > e=rnorm(n) > Z=rep(0,n) > for(t in 3:n) Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+e[t] > Z=Z[800:1000] > n=length(Z) > plot(Z,type="l") Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . Several techniques can be used to estimate those parameters. • using least square regression A natural idea is to see here a regression model, since (if we consider a matrix formulation) Here we can run (conditional) ordinary least squares estimation, > base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)]) > regression=lm(Y~0+X1+X2,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -3.0268 -0.7063 0.1065 0.6925 3.2566 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.23400 0.05463 4.283 2.88e-05 *** X2 0.62863 0.05476 11.479 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.062 on 197 degrees of freedom Multiple R-squared: 0.6349, Adjusted R-squared: 0.6312 F-statistic: 171.3 on 2 and 197 DF, p-value: < 2.2e-16 so we get the following estimators, for the autocorrelation coefficients, and the volatility of the noise > regression$coefficients
X1        X2
0.2339959 0.6286321
> summary(regression)$sigma [1] 1.061839 • using Yule-Walker equations As we’ve seen in class, we can easily get the following equations for the autocovariance functions, which can also be written (again, using a matrix expression) So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions The code is the following > rho1=cor(Z[1:(n-1)],Z[2:n]) > rho2=cor(Z[1:(n-2)],Z[3:n]) > A=matrix(c(1,rho1,rho1,1),2,2) > b=matrix(c(rho1,rho2),2,1) > (PHI=solve(A,b)) [,1] [1,] 0.2256270 [2,] 0.6315329 Now, we need to extract the estimated innovation process, from this set of parameters > estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.058558

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

An alternative could be to include the variance term in Yule-Walker equations, to get a three dimensional linear equation,

$\left\{\begin{array}{l} \gamma_0 = \varphi_1 \gamma_1+\varphi_2 \gamma_2+\sigma^2\\ \gamma_1=\varphi_1 \gamma_0+\varphi_2 \gamma_1 \\ \gamma_2=\varphi_1 \gamma_1+\varphi_2 \gamma_0\end{array}\right.$

It is not much more complicated to solve, actually,

> gamma0=var(Z[1:n])
> gamma1=var(Z[1:(n-1)],Z[2:n])
> gamma2=var(Z[1:(n-2)],Z[3:n])
> A=matrix(c(gamma1,gamma0,gamma1,gamma2,gamma1,gamma0,1,0,0),3,3)
> b=matrix(c(gamma0,gamma1,gamma2),3,1)
> (PHISIGMA=solve(A,b))
[,1]
[1,] 0.2283151
[2,] 0.6283431
[3,] 1.1335501
• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. The standard model is a Gaussian model, i.e.

$Y_t\vert Y_{t-1}=y_{t-1},Y_{t-2}=y_{t-2}$

has a Gaussian distribution

$\mathcal{N}(\varphi_1y_{t-1}+\varphi_2y_{t-2},\sigma^2)$

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par [1] 0.2339589 0.6285002 1.0565613$value
[1] 293.3042

$counts function gradient 106 NA$convergence
[1] 0

$message NULL It is also possible to consider a global maximum likelihood optimisation problem, since the variance matrix of vector $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ has a know form. • using (unconditional) likelihood estimators The variance matrix of $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ is $\boldsymbol{\Gamma}=[\gamma(\vert i-j\vert)]$, where autocovariances are not not know, be can easily be computed using a recursive relationship. > library(mnormt) > GlobalLogLik=function(A,TS){ + n=length(TS) + phi1=A[1]; phi2=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=1 + rho[2]=phi1/(1-phi2) + for(h in 3:n) rho[h]=phi1*rho[h-1]+phi2*rho[h-2] + for(i in 1:n){for(j in 1:n){ + SIG[i,j]=rho[abs(i-j)+1]}} + gamma0=(1-phi2)*sigma^2/((1+phi2)*((1-phi2)^2-phi1^2)) + SIG=gamma0*SIG + return(dmnorm(TS,rep(0,n),SIG,log=TRUE))} > LogL=function(A) -GlobalLogLik(A,TS=Z) > optim(c(.1,.1,1),LogL) Error in chol.default(x, pivot = FALSE) : Error in pd.solve(varcov, log.det = TRUE) : x appears to be not positive definite The problem is that there is a strong constraint on the pair $(\varphi_1,\varphi_2)$ to get a stationary process (we are not far away, here, from the border of the triangle, where the process become non stationary). To be more specific (this was mentioned in a previous post), we should have $\left\{\begin{array}{l} \phi_2-\phi_1<1 \\\phi_2+\phi_1<1\\ \vert\phi_2\vert<1\end{array}\right.$ i.e. in a standard matrix form $\left[\begin{array}{cc} +1 & -1 \\ -1 & -1 \\ 0 & +1\end{array}\right]\left[\begin{array}{c} \varphi_1 \\ \varphi_2\end{array}\right] > \left[\begin{array}{c} -1 \\ -1 \\ -1\end{array}\right]$ (we can add an additional constraint on the variance parameter, to insure that it will be positive). To run a contrained optimization routine, consider > U=matrix(c(1,0,0,-1,0,1,0,-1,0,0,1,0),4,3) > C=c(0,0,0,-.99999) > constrOptim(c(.1,.1,1),LogL,grad=NULL,ui=U,ci=C)$par
[1] 0.2238892 0.6342850 1.0613388

$value [1] 297.9202$counts
108       NA

$convergence [1] 0$message
NULL

$outer.iterations [1] 2$barrier.value
[1] 0.000189892

(here, to faster, we restrain the parameters so that they will be positive).

• comparing those estimates

Here, our five estimators are rather close. Let us run more samples to see more precisely how they behave. For the first parameter $\widehat{\varphi_1}$, we get

and for the second one, $\widehat{\varphi_2}$, we have

The bias we observe is probably coming from the fact that, with this numerical example, we are not far away from the non-stationary case (the sum of the true parameters should be less than 1, and it is 0.95). When we estimate the parameters, we force them to be inside the triangle, since those parameters can be estimated only if the process is stationary.

Observe that the standard-deviation of the innovation process $\widehat{\sigma}$ is here, well estimated,

(with clearly some estimators that perform better than others).

# Bias of Hill Estimators

In the MAT8595 course, we’ve seen yesterday Hill estimator of the tail index. To be more specific, we did see see that if $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, then Hill estimators for $\alpha$ are given by

$\widehat{\alpha}_k = \left[\frac{1}{k}\sum_{i=0}^{k-1} \log X_{n-i,n} -\log X_{n-k,n}\right]^{-1}$
for $k\in\{1,2,\cdots,n\}$. Then we did say that $\widehat{\alpha}_k$ satisfies some consistency in the sense that $\widehat{\alpha}_k \overset{\mathbb{P}}{\rightarrow} \alpha$ if $k\rightarrow\infty$, but not too fast, i.e. $k/n\rightarrow0$ (under additional assumptions on the rate of convergence, it is possible to prove that $\widehat{\alpha}_k \overset{a.s.}{\rightarrow} \alpha$). Further, under additional technical conditions

$\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}(0,\alpha^2)$

In order to illustrate this point, consider the following code. First, let us consider a Pareto survival function, and the associated quantile function

> alpha=1.5
> S=function(x){ifelse(x>1,x^(-alpha),1)}
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root} The code here is obviously too complicated, since this power function can easily be inverted. But later on, we will consider a more complex survival function. Here are the survival function, and the quantile function, > u=seq(0,5,by=.01) > plot(u,Vectorize(S)(u),type="l",col="red") > u=seq(0,99/100,by=.01) > plot(u,Vectorize(Q)(u),type="l",col="blue",ylim=c(0,20)) Here, we need the quantile function to generate a random sample from this distribution, > n=500 > set.seed(1) > X=Vectorize(Q)(runif(n)) Hill plot is here > library(evir) > hill(X) > abline(h=alpha,col="blue") We can now generate thousands of random samples, and see how those estimators behave (for some specific $k$‘s). > ns=10000 > HillK=matrix(NA,ns,10) > for(s in 1:ns){ + X=Vectorize(Q)(runif(n)) + H=hill(X,plot=FALSE) + hillk=function(k) H$y[H$x==k] + HillK[s,]=Vectorize(hillk)(15*(1:10)) + } and if we compute the average, > plot(15*(1:10),apply(HillK,2,mean) we do get a series of estimators that can be considered as unbiased. So far, so good. Now, recall that being in the max-domain of attraction of the Fréchet distribution does not mean that $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, but is means that $\overline{F}(x)= x^{-\alpha} \mathcal{L}(x)$ for some slowly varying function $\mathcal{L}$, not necessarily constant! In order to understand what could happen, we have to be slightly more specific. And this can be done only by looking at second order regular variation property of the survival function. Assume, here that there is some auxilary function $a$ such that $\lim_{t\rightarrow\infty}\frac{\overline{F}(xt)/\overline{F}(t)-x^{-\alpha}}{a(t)}=x^{-\alpha}\frac{1-x^{-\beta}}{\beta}{}$ This (positive) constant $\beta$ is – somehow – related to the speed of convergence of the ratio of the survival functions to the power function (see e.g. Geluk et al. (2000) for some examples). To be more specific, assume that $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ then, the second order regular variation property is obtained using $a(t)=\beta t^{-\beta}$, and then, if $k$ goes to infinity too fast, then the estimator will be biased. More precisely (see Chapter 6 in Embrechts et al. (1997)), if $k=O(n^{2\beta/(\alpha+2\beta)})$, then, for some $\lambda>0$, $\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}\left(\frac{\alpha^3}{\beta-\alpha}\lambda,\alpha^2\right)$ The intuitive interpretation of this result is that if $k$ is too large, and if the underlying distribution is not exactly a Pareto distribution (and we do have this second order property), then Hill’s estimator is biased. This is what we mean when we say • if $k$ is too large, $\widehat{\alpha}_k$ is a biased estimator • if $k$ is too small, $\widehat{\alpha}_k$ is a volatile estimator (the later comes from properties of a sample mean: the more observations, the less the volatility of the mean). Let us run some simulations to get a better understanding of what’s going on. Using the previous code, it is actually extremly simple to generate a random sample with survival function $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ > beta=.5 > S=function(x){+ ifelse(x>1,.5*x^(-alpha)*(1+x^(-beta)),1) } > Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}

If we use the code above. Here, with

> n=500
> set.seed(1)
> X=Vectorize(Q)(runif(n))

the Hill plot becomes

> library(evir)
> hill(X)
> abline(h=alpha,col="blue")

But it’s based on one sample, only. Again, consider thousands of samples, and let us see how Hill’s estimator is behaving,

so that the (empirical) mean of those estimator is

# Informatique (sans ordinateur), partie 1

Pendant les vacances, avec mon fils, on s’est amusé à faire les premières activités proposés dans le manuel Informatique Sans Ordinateur, disponible sur http://csunplugged.org/ (où plusieurs activités complémentaires peuvent être téléchargées, mais en anglais seulement). Pour des raisons personnelles, mon fils est de plus en plus amené à utiliser l’ordinateur à l’école. Et je dois avouer que ça me gêne de la voir manipuler un outils qu’il ne maîtrise pas vraiment (je renvoie d’ailleurs à un billet de Dr Goulu qui disait la même chose au début de l’année). Que les choses soient claires : je ne prétend pas non plus maîtriser l’informatique ! Mais comme je l’ai déjà raconté sur ce blog, à son âge, je codais mes premiers jeux (souvent en recopiant des codes en BASIC), et j’ai l’impression que manipulant les ordinateurs depuis presque 30 ans m’a permis d’avoir un peu de recul, en tous les cas plus que lui. Et je dois avouer que le manque de culture en informatique de la génération entre la mienne et celle de mon fils me surprend. Bref, comme je l’avais dit dans un précédant billet, ça me dérange qu’une génération autant amenée à utiliser l’outil informatique en sache aussi peu à ce sujet. Aussi, quand j’ai découvert ce petit programme d’activités le mois dernier, j’ai eu envie de tenter l’expérience.

Au départ, je pensais proposer les activités à mon fils et à ma fille (qui a 8 ans) mais cette dernière était plongée dans des dessins quand on a débuté, et je n’ai pas réussi à la faire décrocher. Pour ceux qui veulent tenter l’expérience, il y a une logique dans les activités, et je pense qu’il serait dommage se rater les premières. Bref, je n’ai pas re-proposé à ma fille de se joindre à nous. Mais on verra aux prochaines vacances….

Le document propose une douzaine d’activités (davantage si on traîne sur http://csunplugged.org/activities). La première partie (que je vais évoquer ici) porte sur la représentation de l’information.

• section 1 : le système binaire

La première activité est vraiment bien faite. On apprend la base 2, l’écriture en 0 et en 1, les codes ASCII, et les notions de 32 et 64 bits. Histoire d’illustrer, je vais coder un peu (sur ordinateur cette fois) pour expliquer comment ça fonctionne.

> base2=function(x,n=8){
+ Base.b=rep(0,n)
+ ndigits=(floor(logb(x,base=2))+1)
+ for(i in 1:ndigits){
+ Base.b[n-i+1]=(x%%2)
+ x=(x %/% 2)}
+ plot(0:1,0:1,xlab="",ylab="",
+ axes=FALSE, xlim=c(0,n),ylim=c(0,1),col="white")
+ for(i in 1:n){
+ polygon(i-1+c(.1,.1,.9,.9),c(.1,.9,.9,.1),lwd=2,
+ col=c("white","red")[1+(Base.b[i]==1)])}
+ return(Base.b)}

On découpe des bouts de cartons, et on les juxtapose, de manière à écrire des nombres.

Regardons par exemple le nombre 17. Rouge pour un 1, et blanc pour un 0.

> base2(17)
[1] 0 0 0 1 0 0 0 1

Le premier jeu est d’écrire quelques nombres, pour se familiariser. Ensuite, on voit que c’est même facile de faire des opérations de base. Par exemple, si on multiplie par 2, c’est facile : on décale vers la gauche, et on rajoute un carré blanc, tout à droite

> base2(17*2)
[1] 0 0 1 0 0 0 1 0

Amusant, non ? On peut aussi faire des additions, par exemple, 12 s’écrit

> base2(12)
[1] 0 0 0 0 1 1 0 0

et si on somme 12 et 17, on obtient

> base2(12+17)
[1] 0 0 0 1 1 1 0 1

où on raisonne comme en base 10. En fait, dans cet exemple, on retrouve que 0+1=1+0=1 et que 0+0=0 (oui, il n’y a pas de retenue dans cet exemple). Ensuite, on utilise une décomposition en nombre de l’alphabet pour coder des lettres (A-1, B-2, C-3, etc), puis rapidement, on passe aux codes ASCII, C’est très ludique ! Je dois avouer qu’on en a profité pour faire une digression par les codes secrets, mais j’en reparlerais une autre fois. Bref, cette première activité nous a emballé !

• section 2 : pixeliser et dessiner.

Suite logique de l’étape précédente, on a ensuite évoqué la pixellisation. Cette activité me parlait beaucoup, car c’est précisément ce que je faisais quand j”avais l’âge de mon fils (ou presque). Sur le MO5 de la maison, il y a avait essentiellement deux jeux, sur cassettes: un jeu de voitures, et super tennis. Au tennis, les personnages étaient assez simples avec 4 positions possibles: l’attente, le service, le coup droit, et le revers. Ensuite, la figure était juste translatée (on avait déjà évoqué ce point il y a quelques mois).

Mais autant ça me parle énormément, autant je me suis demandé ce que ça pouvait dire aujourd’hui, car les images sont incroyablement lissées… On a pris une image sur mon ordinateur, et on a zoomé, zoomé…

… en vain. On ne voit plus vraiment les pixels sur les dessins. C’est dommage ! Heureusement, j’ai fini par trouver quelques exemple pour illustrer ce concept (finalement assez théorique) de la pixellisation. Je ne me sentais pas très à l’aise pour parler de lissage avec mon fils. Ce sont des choses que j’explique à l’occasion à mes étudiants, et ça me gêne que mon fils de 11 ans sache plus de choses que mes étudiants. En fait, ce dilemme m’a tiraillé tout au long de nos activités, pour être honnête… D’un côté, ce sont des exercices pour enfants de primaire, et de l’autre, ce sont des choses que certains étudiants en formation universitaire de mathématique se devraient de maîtriser.

Pour revenir un peu à l’activité, c’était amusant, et ça annonçait un peu le principe de la compression (qui viendrait avec la prochaine activité). Mais trouver une tasse, et une image de Saturne en coloriant des cases, ça a vite lassé mon fils. On a fini par abandonner sans vraiment finir l’activité. Cela dit, j’ai découvert par la suite qu’il y avait des activités amusantes en ligne, sur http://csunplugged.org/activities, avec en particulier une discussion sur les traits et les cercles (en anglais seulement)

• section 3 : compresser et zipper

Cette activité était…. surprenante. On a joué avec des algorithmes de compression, type LZ77 utilisé pour zipper des fichiers. Le principe est assez amusant…. Par exemple. dans la phrase

on retrouve des blocs de lettres. Alors l’idée (difficile à suivre si on se limite aux indications données) est que si on a plus de 2 lettres qui se répètent, on peut faire un pointage.

On va alors remplacer le second bloc par de l’information expliquant aller chercher l’information,

Aussi, la phrase devient

Autrement dit, on lui dit de reculer d’un certain nombre de caractères, et d’utiliser un bloc de caractères, d’une certaine longueur. Amusant non ? En plus, on peut utiliser une écriture récursive, et remplacer par des lettres non encore définies. Comme le joli exemple suivant

On recule de 2 lettres, et on prend un bloc de 3 lettres. Amusant non ? C’est facile à faire, et on peut facilement lire un fichier compressé. Par contre, compresser soit même un texte est plus compliqué ! On s’est longtemps interrogé sur la marche à suivre, car l’aogorithme est mal expliqué. Manifestement, on part de la gauche, et si on retrouve un bloc de 2 caractères déjà lu (ou plus), on pointe dessus. Manifestement, selon ce qu’on y lit, on peut passer de 2500 caractères à 500. C’est dommage que la compression en question ne soit pas mentionnée. Par contre, on imagine bien que plus le texte est long, plus on va retrouver des blocs déjà vus. Et l’idée d’illustrer par un poème est brillante !

On s’est bien amusé pendant cette activité, on a fait des essais, mais on ne sait pas si ce qu’on fait est “optimal“. Cela dit,  m’a donné envie de lire davantage, d’autant que plusieurs articles évoquent ce point (et que la notion d’information et d’entropie se cachent derrière, mais on va la voir dans quelques sections).

• section 4 : gérer les erreurs et l’exemple des codes barres

On a ensuite une activité réellement amusante, sur la gestion des erreurs. Ça commence par un petit tour de magie et l’écriture binaire…. On considère la grille suivante (avec les cartes utilisées dans la première activité: une face rouge, une face blanche)

Les couleurs ont étés mises au hasard… on a constitué une grille 5×5 en posant les cartes au hasard (ce sont les cartes sur le tapis bleu). Pour ceux qui sont autour… c’est le truc du tour de magie. Maintenant, je lasse mon fils retourner une carte au hasard, sur le tapis bleu (en fait, je pense qu’il peut retourner celle qu’il veut)

Pendant qu’il la retournait, je n’ai pas regardé, et il va falloir maintenant retrouver quelle carte a été retournée… En fait, le truc pour la retrouver, c’est que les cartes à droite, et en bas, ont été mises de telle sorte que le nombre de cartes rouges par ligne, et par colonne soit pair !

En plus, si deux cartes avaient été retournées, j’aurais aussi pu les retrouver. Amusant non ? L’idée est vraiment géniale…. Ensuite, on apprend à vérifier les codes d’identité bancaire ou de code ISBN, pour les livres. Oui, le code qui figure sur les livres, au dos.

Bon, le soucis est que sur ce code barre, il y a manifestement un problème, car le code généré n’est pas le bon. En fait, le code ISBN qu’on va utiliser est le suivant

> isbn=1466592591
> checkcode=function(n){
+ l=as.character(n)
+ while(nchar(l)<10) l=paste("0",l,sep="")
+ a=substr(l,1,9)
+ y=as.numeric(substr(l,10,10))
+ x=as.numeric(unlist(strsplit(a,"")))
+ s=sum(x*(10:2))
+ z=11-(s%%11)
+ return(z==y)
+ }

Sur notre code à 10 chiffres, le dernier chiffre va servir à vérifier si les premiers sont bons, ou pas. Tout est expliqué dans la définition du International Standard Book Number sur Wikipedia. C’est presque simple, à condition de maîtriser le reste de la division par 11 (comme le montre le code précédant)

> checkcode(isbn)
[1] TRUE

On a pris quelques livres dans la bibliothèque, et je lui demandais de me dicter le code…. C’était assez amusant comme exemple… Mais mon fils a voulu essayer un nombre au hasard… et il a trouvé un code ISBN du premier coup…. Bon, il avait une chance sur 12.

> mean(Vectorize(checkcode)(trunc(1e10*runif(1000000))))
[1] 0.08327
• section 5 : jeu du “devine le chiffre auquel je pense” et système binaire

La dernière section a été de loin la plus intéressante ! On effleure l’information (au sens de Shannon) et le logarithme (en base 2)… J’avais évoqué des idées similaires dans un vieux billet, justement, suite à un jeu avec les enfants. On voit en particulier la construction par arbre de la méthode dichotomique… Par exemple, si je pense à un nombre (entier) entre 0 et 7, la méthode pour trouver le nombre, la plus efficace, est la suivante :

Ce qui est amusant, c’est qu’on retrouve la décomposition en base binaire des nombres,

avec d’abord 2² puis 2¹ et enfin 2⁰. Amusant, non ? On a ensuite fait un test (et ma fille est venu se joindre à nous). On a commencé par “choisir un nombre entre 1 et 100” puis “choisir un nombre entre 1 et 1000“, à tour de rôle. On a vu qu’il fallait environ 7 coups pour réussir avec 100 chiffres possibles, et une dizaine avec 1000. Sans parler de logarithme, on a vu qu’il fallait chercher les exposants de 2 pour attendre 100, ou 1000.

Mais le plus intéressant, c’est que j’ai commencé, pour montrer qu’on divise en deux, à chaque fois. Je commençais avec “plus grand que 50 ?” puis “plus grand que 25 ?” puis “plus grand que 12 ?” etc. Mon fils a opté pour une stratégie plus étrange, en commençant par couper au milieu “plus grand que 50 ?” puis (directement) “plus grand que 5 ?“. En fait, en jouant plusieurs fois, je me suis rendu compte que sa stratégie pouvait être (en quelque sorte) optimale, car sa sœur ne prend visiblement pas des nombres de manière uniforme entre les deux bornes: elle a tendance à prendre des petits ou des très grands nombres. Sur une partie, il a ainsi être plus rapide que ma stratégie soit-disant optimale, en gagnant en 5 coups. On a ensuite eu une longue discussion sur ce que pouvait être la stratégie optimale. C’était vraiment une discussion intéressante… qui s’est poursuivie avec l’activité suivante, de classement et de tri. Mais on en reparlera bientôt !

Quand on voit les ressources qui traînent en ligne, je ne comprends pas que l’enseignement de l’informatique en primaire ne soit pas obligatoire : c’est incroyablement ludique ! Bref, on s’est vraiment bien amusé avec ces activités, et on a appris plein de choses. Bon, peut être que le fait d’être resté autour de -30°C ces derniers jours n’a pas vraiment permis de proposer d’alternatives…

# Modèle de régression et interaction(s) entre facteurs

Dans un modèle de régression, on veut écrire

$\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\varphi(\boldsymbol{X})=\varphi(X_1,\ldots,X_d)$

Quand on se limite à un modèle linéaire, on écrit

$\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta}$

ou encore

$\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\beta_0+\sum_{i=1}^d \beta_i X_i$

Mais on de doute que l’on rate quelque chose… en particulier, on va rater toutes les interactions possibles. On peut croiser les variables, et supposer que

$\mathbb{E}(Y\vert \boldsymbol{X}=(X_1,\ldots,X_d))=\beta_0+\sum_{i=1}^d \beta_i X_i+\sum_{i=1}^{d-1}\sum_{j=i+1}^d \gamma_{i,j} X_i X_j$

qui peut s’étendre d’avantage, à l’ordre 3,

$\beta_0+\sum_{i=1}^d \beta_i X_i+\sum_{i=1}^{d-1}\sum_{j=i+1}^d \gamma_{i,j} X_i X_j+\sum_{i=1}^{d-2}\sum_{j=i+1}^{d-1}\sum_{k=j+1}^{d} \delta_{i,j,k} X_i X_j X_k$

voire davantage.

Supposons que nos variables $X_i$ soient ici qualitatives, et plus précisément binaires. Prenons un exemple simple, avec des données (classiques) en risque de crédit1. On peut trouver la base via

library(evtree)
db=GermanCredit

ou encore directement

myVariableNames = c("checking_status","duration","credit_history",
"purpose","credit_amount","savings","employment","installment_rate",
"personal_status","other_parties","residence_since","property_magnitude",
"age","other_payment_plans","housing","existing_credits","job",
"num_dependents","telephone","foreign_worker","class")

"http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",
header=FALSE,col.names=myVariableNames)

Retenons pour commencer trois variables explicatives,

db=data.frame(Y=GermanCredit$class-1, X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"), X3=GermanCredit$savings%in%c("A61","A62"))
reg=glm(Y~X1+X2+X3,data=db,family=binomial)
summary(reg)

La régression sans interaction donne ici

Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.5431  -0.8421  -0.6295   1.3994   1.9999

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.8544     0.1699 -10.915  < 2e-16 ***
X1TRUE        0.3363     0.1496   2.249   0.0245 *
X2TRUE        1.3462     0.2347   5.735 9.76e-09 ***
X3TRUE        1.0001     0.1787   5.596 2.19e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1221.7  on 999  degrees of freedom
Residual deviance: 1143.6  on 996  degrees of freedom
AIC: 1151.6

Number of Fisher Scoring iterations: 4

Il existe plusieurs interactions possibles ici (limitons nous aux paires). C’est ce que l’on observe quand on fait la régression

reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3, family = binomial,
data = db)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.5369  -0.8281  -0.6439   1.3954   1.9638

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.77109    0.20070  -8.825  < 2e-16 ***
X1TRUE         0.30296    0.33737   0.898 0.369186
X2TRUE         0.88353    0.54255   1.628 0.103421
X3TRUE         0.87709    0.22583   3.884 0.000103 ***
X1TRUE:X2TRUE -0.37917    0.49343  -0.768 0.442225
X1TRUE:X3TRUE  0.09178    0.37278   0.246 0.805522
X2TRUE:X3TRUE  0.80923    0.58185   1.391 0.164293
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1221.7  on 999  degrees of freedom
Residual deviance: 1141.0  on 993  degrees of freedom
AIC: 1155

Number of Fisher Scoring iterations: 4

On peut faire un dessin pour visualiser les interactions : on a trois sommets (nos trois variables), et on visualiser les interactions

indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",
xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

ce qui donne ici, pour nos trois variables

Ce modèle pourrait sembler incomplet, car on ne regarde que les interactions entre les modalités, par paires. En fait, c’est parce qu’il manque (visuellement) les variables non-croisées. On peut les rajouter si on veut (au risque de surcharger le dessin)

cercle=function(c,r,cl) lines(c[1]+r*cos(seq(0,2*pi,length=501)),
c[2]+r*sin(seq(0,2*pi,length=501)),col=cl)

reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

soit ici

Si on change le ‘sens‘ de nos variables (en recodant a l’envers, en permutant les vrais et les faux), on obtient le graphique suivant

dbinv=db
dbinv[,2:k]=1-dbinv[,2:k]
reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=dbinv,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

qui peut alors être comparé au graphique précédant

Avec 5 variables, on augmente les interactions possibles… même si beaucoup risquent d’être non-significatifs. On peut déjà se focaliser sur les paires possibles d’interactions croisées. Pour simplifier le code, on va utiliser deux fonctions locales,

vrepeach=function(x,e){
v=NULL
for(i in 1:length(e)){v=c(v,rep(x[i],each=e[i]))}
return(v)}
vreplength=function(x,l){
v=NULL
for(i in 1:length(l)){v=c(v,x[l[i]:length(x)])}
return(v)}

et ensuite, on adapte le code précédant

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

ce qui donne un schéma plus complexe,

On peut aussi prendre juste 2 variables, prenant 3 et 4 modalités respectivement. On va extraire deux variables indicatrices pour la première (la modalité restante sera la modalité de référence) et trois pour la seconde,

db=data.frame(Y=GermanCredit$class-1, X1=GermanCredit$checking_status=="A12",
X2=GermanCredit$checking_status=="A13", X3=GermanCredit$checking_status=="A14",
X4=GermanCredit$employment%in%c("A72","A73"), X5=GermanCredit$employment%in%c("A74","A75"))
k=5
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

On voit que plusieurs interactions ne sont alors plus possibles, sur la partie gauche (les trois modalités de la même variable) et sur la partie droite

On peut d’ailleurs simplifier les graphs, en ne visualisant que les interactions significatives.

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) }}} for(i in 1:k){ if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

soit ici

Ici, une seule interactions croisée est significative, et presque toutes les variables le sont. Et si on reprend le modèle avec 5 facteurs,

db=data.frame(Y=GermanCredit$class-1,X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"), X3=GermanCredit$savings%in%c("A61","A62"),
X4=GermanCredit$employment%in%c("A71","A72"), X5=GermanCredit$other_payment_plans=="A143")

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){ segments(sommetX[indices[i,2]],sommetY[indices[i,2]], sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey") text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]), trunc(10000*coefficients(reg)[1+k+i])/10000) }}} for(i in 1:k){ if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

on obtient

Je ne sais pas si mes graphiques sont pertinents, ou pas. Mais je trouve ça joli. En fait, je suis tombé un peu par hasard2 sur les Tables de Taguchi, développées par Gen’ichi Taguchi (田口 玄一). Le soucis est que je n’ai rien compris… Enfin, disons que je croyais comprendre, puis j’ai continué à faire des dessins… Si quelqu’un pourrait m’expliquer sur mon exemple les graphiques de Taguchi, je suis preneur ! car je doute que ce soit ce que je fais depuis tout à l’heure…

1. Cette base est largement utilisée dans le quatrième chapitre de Computational Actuarial Science with R, à paraître dans les mois à venir.

2.En l’occurence, le hasard est @Benavent qui a suscité ma curiosité ce matin en me parlant de ces tables, dont je n’avais alors jamais entendu parlé ! J’avais même lu rapidement Taniguchi (谷口 ジロー) et je ne voyais pas le rapport avec les statistiques….

# Pricing reinsurance contracts, another case study

A reinsurance case study for tomorrow’s class. The goal will be to price some nonproportional reinsurance contract, for business interruption claims. Consider the following dataset,

> library(gdata)
+ "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls",
+  sheet=1)
Content type 'application/vnd.ms-excel' length 183808 bytes (179 Kb)
open URL
==================================================
downloaded 179 Kb

As for any (standard) insurance contract, there are two parts in the pricing

• the expected number of claims
• the average cost of individual claims

Here, we do not have covariates (but it might be possible to use some, like the kind of industry, the location, etc).

Let us start with the expected number of claims, per year. Here is the daily frequency,

The data are rather old… but somehow, it is a good thing since after ten years, we can expect that most of the claims have been settled (we’ll discuss claims dynamic starting next week). To plot the graph above, we use

> date=db$DSUR > D=as.Date(as.character(date),format="%Y%m%d") > vD=seq(min(D),max(D),by=1) > sD=table(D) > d1=as.Date(names(sD)) > d2=vD[-which(vD%in%d1)] > vecteur.date=c(d1,d2) > vecteur.cpte=c(as.numeric(sD),rep(0,length(d2))) > base=data.frame(date=vecteur.date,cpte=vecteur.cpte) > plot(vecteur.date,vecteur.cpte,type="h",xlim=as.Date(as.character( + c(19850101,20111231)),format="%Y%m%d")) Then, we can get a prediction of the daily number of business interruption claims, e.g. for any day in 2010 (assume that we had to price a reinsurance contract a few years ago), using a (standard) Poisson regression > regdate=glm(cpte~date,data=base,family=poisson(link="log")) > nd2010=data.frame(date=seq(as.Date(as.character(20100101),format="%Y%m%d"), + as.Date(as.character(20101231),format="%Y%m%d"),by=1)) > pred2010 =predict(regdate,newdata=nd2010,type="response") > sum(pred2010) [1] 159.4757 Observe that using old data has drawbacks, since we got much more uncertainty if we use a regression on time (to include some possible trend) Say we have something like 160 claims over a given year, on average. > plot(D,db$COUTSIN,type="h")

Let us now focus on the cost of those claims. We have 2,400 claims in our dataset, to fit a model (or at least estimate how much a reinsurance contract might cost us). Assume that we would like to purchase a reinsurance contract for our very large claims. Like the two largest per year. Over 16 years, the decutible should be close to the cost of the 32nd largest claim, which was close to 15 million.

> quantile(db$COUTSIN,1-32/2400)/1e6 98.66667% 15.34579 > abline(h=quantile(db$COUTSIN,1-32/2400),col="blue")

So consider some reinsurance contract with a deductible of 15 million. Unfortunately, we cannot find unlimited covers. So let us assume that a reinsurance company agrees for such a deductible, but with a limited cover of 35 million. The average cost (for the reinsurance company) is $\mathbb{E}(g(X))$ where

$g(x)=\min\{35,\max\{x-15,0\}\}$

A first idea is to look at the first cost, i.e. the empirical average of that indemnity, on our portfolio. The indemnity function is

> indemn=function(x) pmin((x-15)*(x>15),50-15)

we can check on a few losses that it is actually what we wish to compute

> indemn(5)
[1] 0
> indemn(20)
[1] 5
> indemn(50)
[1] 35

Now, if the compute the average repayment by the reinsurance company, over 16 years, we get

> mean(indemn(db$COUTSIN/1e6)) [1] 0.1624292 So, per claim, the reinsurance company will pay, on average 162,430. With 160 claims per year, the pure premium should be close to 26 million > mean(indemn(db$COUTSIN/1e6))*160
[1] 25.98867

(again, for a 35 million cover, for some claims that should occur, on average, twice a year). As we will see, a standard model in reinsurance is the Pareto distribution (or to be more specific, a Generalized Pareto one),

$F_{(\xi,\mu,\sigma)}(x) = \begin{cases}1 - \left(\ds{1+ \frac{\xi(x-\mu)}{\sigma}\right)^{-1/\xi}} & \text{for }\xi \neq 0, \\1 - \exp \left(-\frac{x-\mu}{\sigma}\right) & \text{for }\xi = 0.\end{cases}$

There are three parameters here

• the threshold $\mu$ (that we will consider as fixed, but we will see its impact on reinsurance pricing)
• the scale parameter $\sigma$ (called $\beta$ in R)
• the tail index $\xi$

The strategy is to consider a threshold below our deductible, e.g. 12 million. Then, given that the loss exceed 12 million, we can fit a Generalized Pareto distribution,

> gpd.PL <- gpd(db$COUTSIN,12e6)$par.ests
> gpd.PL
xi         beta
7.004147e-01 4.400115e+06

and compute

$\int_{d}^{d+c}(x-d)\ dF_{(\xi,\mu,\sigma)}(x)$

>  E <- function(yinf,ysup,xi,beta,threshold){
+    as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta),
+    lower=yinf,upper=ysup)$value+ + (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf)) + } Here, given that a claim exceeds 12 million, the average repayment is close to 6 million > E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 6058125 Now, we have to take into account the probability to reach 12 million, which is here > mean(db$COUTSIN>12e6)
[1] 0.02639296

So, if we summarize, we have on average 160 claims per year,

> p
[1] 159.4757

Only 2.6% will exceed 12 million

> mean(db$COUTSIN>12e6) [1] 0.02639296 So, the yearly frequency of claism larger than 12 million is 4.2 claims > p*mean(db$COUTSIN>12e6)
[1] 4.209036

And for a claim that exceed 12 million, the average repayment is

> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6)
[1] 6058125

So, the pure premium should be close to

> p*mean(db$COUTSIN>12e6)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 25498867 which (hopefully) is close to the empirical value we got. Actually, it is also possible to look at the impact of the threshold parameter, since it is clearly and intermediate value that could be changed. I mean, why 12 and not 10? Consider > esp=function(threshold=12e6,p=sum(pred2010)){ + (gpd.PL <- gpd(db$COUTSIN,threshold)$par.ests) + return(p*mean(db$COUTSIN>threshold)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],threshold))
+  }

We can plot the pure premium as a function of that threshold,

> seuils=seq(1e6,15e6,by=1e6)
> plot(seuils,Vectorize(esp)(seuils),type="b",col="red")

which is between 24 and 26 for large thresholds. Again, that is only the first step, and we can price a higher reinsurance layer, like a reinsurance contract with a deductible of 50 million (we have our previous reinsurance contract for claims below that threshold), and a cover of 50 million, for instance. For those high layers, it become interesting to have a parametric model, which should be more robust than the empirical average.

# Generating your own normal distribution table

It might sounds incredibly old fashion, but for my the exam for the ACT2121 probability course (to prepare for the exam P of the Society of Actuaries), I will provide a standard normal distribution table. The problem is that it is never the one we’re looking for (sometimes it is the survival function, sometimes it is the cumulative distribution function, sometimes we consider only positive values, etc). Here is the one that will be given for the exam, this Friday.

Now, here is the code to generate it.

I did use the following code to generate the table (in a latex format),

> u=seq(0,3.09,by=0.01)
> p=pnorm(u)
> m=matrix(p,ncol=10,byrow=TRUE

We have here the table that we wish to have in our table,

> options(digits=4)
> m
[,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
[1,] 0.5000 0.5040 0.5080 0.5120 0.5160 0.5199 0.5239 0.5279 0.5319 0.5359
[2,] 0.5398 0.5438 0.5478 0.5517 0.5557 0.5596 0.5636 0.5675 0.5714 0.5753
[3,] 0.5793 0.5832 0.5871 0.5910 0.5948 0.5987 0.6026 0.6064 0.6103 0.6141
[4,] 0.6179 0.6217 0.6255 0.6293 0.6331 0.6368 0.6406 0.6443 0.6480 0.6517
[5,] 0.6554 0.6591 0.6628 0.6664 0.6700 0.6736 0.6772 0.6808 0.6844 0.6879
[6,] 0.6915 0.6950 0.6985 0.7019 0.7054 0.7088 0.7123 0.7157 0.7190 0.7224
[7,] 0.7257 0.7291 0.7324 0.7357 0.7389 0.7422 0.7454 0.7486 0.7517 0.7549
[8,] 0.7580 0.7611 0.7642 0.7673 0.7704 0.7734 0.7764 0.7794 0.7823 0.7852
[9,] 0.7881 0.7910 0.7939 0.7967 0.7995 0.8023 0.8051 0.8078 0.8106 0.8133
[10,] 0.8159 0.8186 0.8212 0.8238 0.8264 0.8289 0.8315 0.8340 0.8365 0.8389
[11,] 0.8413 0.8438 0.8461 0.8485 0.8508 0.8531 0.8554 0.8577 0.8599 0.8621
[12,] 0.8643 0.8665 0.8686 0.8708 0.8729 0.8749 0.8770 0.8790 0.8810 0.8830
[13,] 0.8849 0.8869 0.8888 0.8907 0.8925 0.8944 0.8962 0.8980 0.8997 0.9015
[14,] 0.9032 0.9049 0.9066 0.9082 0.9099 0.9115 0.9131 0.9147 0.9162 0.9177
[15,] 0.9192 0.9207 0.9222 0.9236 0.9251 0.9265 0.9279 0.9292 0.9306 0.9319
[16,] 0.9332 0.9345 0.9357 0.9370 0.9382 0.9394 0.9406 0.9418 0.9429 0.9441
[17,] 0.9452 0.9463 0.9474 0.9484 0.9495 0.9505 0.9515 0.9525 0.9535 0.9545
[18,] 0.9554 0.9564 0.9573 0.9582 0.9591 0.9599 0.9608 0.9616 0.9625 0.9633
[19,] 0.9641 0.9649 0.9656 0.9664 0.9671 0.9678 0.9686 0.9693 0.9699 0.9706
[20,] 0.9713 0.9719 0.9726 0.9732 0.9738 0.9744 0.9750 0.9756 0.9761 0.9767
[21,] 0.9772 0.9778 0.9783 0.9788 0.9793 0.9798 0.9803 0.9808 0.9812 0.9817
[22,] 0.9821 0.9826 0.9830 0.9834 0.9838 0.9842 0.9846 0.9850 0.9854 0.9857
[23,] 0.9861 0.9864 0.9868 0.9871 0.9875 0.9878 0.9881 0.9884 0.9887 0.9890
[24,] 0.9893 0.9896 0.9898 0.9901 0.9904 0.9906 0.9909 0.9911 0.9913 0.9916
[25,] 0.9918 0.9920 0.9922 0.9925 0.9927 0.9929 0.9931 0.9932 0.9934 0.9936
[26,] 0.9938 0.9940 0.9941 0.9943 0.9945 0.9946 0.9948 0.9949 0.9951 0.9952
[27,] 0.9953 0.9955 0.9956 0.9957 0.9959 0.9960 0.9961 0.9962 0.9963 0.9964
[28,] 0.9965 0.9966 0.9967 0.9968 0.9969 0.9970 0.9971 0.9972 0.9973 0.9974
[29,] 0.9974 0.9975 0.9976 0.9977 0.9977 0.9978 0.9979 0.9979 0.9980 0.9981
[30,] 0.9981 0.9982 0.9982 0.9983 0.9984 0.9984 0.9985 0.9985 0.9986 0.9986
[31,] 0.9987 0.9987 0.9987 0.9988 0.9988 0.9989 0.9989 0.9989 0.9990 0.9990
> rownames(m)=seq(0,3,b=.1)
> colnames(m)=seq(0,.09,by=.01)

To put it in a nice latex format, we can use

> library(xtable)
> newm=xtable(m,digits=4)
> print.xtable(newm, type="latex", file="nor1.tex")

We now have a simple tex file containing a table.

\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrrrrrr}
\hline
& 0 & 0.001 & 0.002 & 0.003 & 0.004 & 0.005 & 0.006 & 0.007 & 0.008 & 0.009 \\
\hline
0 & 0.5000 & 0.5040 & 0.5080 & 0.5120 & 0.5160 & 0.5199 & 0.5239 & 0.5279 & 0.5319 & 0.5359 \\
0.1 & 0.5398 & 0.5438 & 0.5478 & 0.5517 & 0.5557 & 0.5596 & 0.5636 & 0.5675 & 0.5714 & 0.5753 \\
0.2 & 0.5793 & 0.5832 & 0.5871 & 0.5910 & 0.5948 & 0.5987 & 0.6026 & 0.6064 & 0.6103 & 0.6141 \\
0.3 & 0.6179 & 0.6217 & 0.6255 & 0.6293 & 0.6331 & 0.6368 & 0.6406 & 0.6443 & 0.6480 & 0.6517 \\
0.4 & 0.6554 & 0.6591 & 0.6628 & 0.6664 & 0.6700 & 0.6736 & 0.6772 & 0.6808 & 0.6844 & 0.6879 \\
0.5 & 0.6915 & 0.6950 & 0.6985 & 0.7019 & 0.7054 & 0.7088 & 0.7123 & 0.7157 & 0.7190 & 0.7224 \\
0.6 & 0.7257 & 0.7291 & 0.7324 & 0.7357 & 0.7389 & 0.7422 & 0.7454 & 0.7486 & 0.7517 & 0.7549 \\
0.7 & 0.7580 & 0.7611 & 0.7642 & 0.7673 & 0.7704 & 0.7734 & 0.7764 & 0.7794 & 0.7823 & 0.7852 \\
0.8 & 0.7881 & 0.7910 & 0.7939 & 0.7967 & 0.7995 & 0.8023 & 0.8051 & 0.8078 & 0.8106 & 0.8133 \\
0.9 & 0.8159 & 0.8186 & 0.8212 & 0.8238 & 0.8264 & 0.8289 & 0.8315 & 0.8340 & 0.8365 & 0.8389 \\
1 & 0.8413 & 0.8438 & 0.8461 & 0.8485 & 0.8508 & 0.8531 & 0.8554 & 0.8577 & 0.8599 & 0.8621 \\
1.1 & 0.8643 & 0.8665 & 0.8686 & 0.8708 & 0.8729 & 0.8749 & 0.8770 & 0.8790 & 0.8810 & 0.8830 \\
1.2 & 0.8849 & 0.8869 & 0.8888 & 0.8907 & 0.8925 & 0.8944 & 0.8962 & 0.8980 & 0.8997 & 0.9015 \\
1.3 & 0.9032 & 0.9049 & 0.9066 & 0.9082 & 0.9099 & 0.9115 & 0.9131 & 0.9147 & 0.9162 & 0.9177 \\
1.4 & 0.9192 & 0.9207 & 0.9222 & 0.9236 & 0.9251 & 0.9265 & 0.9279 & 0.9292 & 0.9306 & 0.9319 \\
1.5 & 0.9332 & 0.9345 & 0.9357 & 0.9370 & 0.9382 & 0.9394 & 0.9406 & 0.9418 & 0.9429 & 0.9441 \\
1.6 & 0.9452 & 0.9463 & 0.9474 & 0.9484 & 0.9495 & 0.9505 & 0.9515 & 0.9525 & 0.9535 & 0.9545 \\
1.7 & 0.9554 & 0.9564 & 0.9573 & 0.9582 & 0.9591 & 0.9599 & 0.9608 & 0.9616 & 0.9625 & 0.9633 \\
1.8 & 0.9641 & 0.9649 & 0.9656 & 0.9664 & 0.9671 & 0.9678 & 0.9686 & 0.9693 & 0.9699 & 0.9706 \\
1.9 & 0.9713 & 0.9719 & 0.9726 & 0.9732 & 0.9738 & 0.9744 & 0.9750 & 0.9756 & 0.9761 & 0.9767 \\
2 & 0.9772 & 0.9778 & 0.9783 & 0.9788 & 0.9793 & 0.9798 & 0.9803 & 0.9808 & 0.9812 & 0.9817 \\
2.1 & 0.9821 & 0.9826 & 0.9830 & 0.9834 & 0.9838 & 0.9842 & 0.9846 & 0.9850 & 0.9854 & 0.9857 \\
2.2 & 0.9861 & 0.9864 & 0.9868 & 0.9871 & 0.9875 & 0.9878 & 0.9881 & 0.9884 & 0.9887 & 0.9890 \\
2.3 & 0.9893 & 0.9896 & 0.9898 & 0.9901 & 0.9904 & 0.9906 & 0.9909 & 0.9911 & 0.9913 & 0.9916 \\
2.4 & 0.9918 & 0.9920 & 0.9922 & 0.9925 & 0.9927 & 0.9929 & 0.9931 & 0.9932 & 0.9934 & 0.9936 \\
2.5 & 0.9938 & 0.9940 & 0.9941 & 0.9943 & 0.9945 & 0.9946 & 0.9948 & 0.9949 & 0.9951 & 0.9952 \\
2.6 & 0.9953 & 0.9955 & 0.9956 & 0.9957 & 0.9959 & 0.9960 & 0.9961 & 0.9962 & 0.9963 & 0.9964 \\
2.7 & 0.9965 & 0.9966 & 0.9967 & 0.9968 & 0.9969 & 0.9970 & 0.9971 & 0.9972 & 0.9973 & 0.9974 \\
2.8 & 0.9974 & 0.9975 & 0.9976 & 0.9977 & 0.9977 & 0.9978 & 0.9979 & 0.9979 & 0.9980 & 0.9981 \\
2.9 & 0.9981 & 0.9982 & 0.9982 & 0.9983 & 0.9984 & 0.9984 & 0.9985 & 0.9985 & 0.9986 & 0.9986 \\
3 & 0.9987 & 0.9987 & 0.9987 & 0.9988 & 0.9988 & 0.9989 & 0.9989 & 0.9989 & 0.9990 & 0.9990 \\
\hline
\end{tabular}
\end{table}

and the following code to get a graph, illustrating was was actually computed, in the table (see a previous post for more details)

> library("tikzDevice")
> options(tikzMetricPackages = c("\\usepackage[utf8]{inputenc}",
+ "\\usepackage[T1]{fontenc}", "\\usetikzlibrary{calc}", "\\usepackage{amssymb}"))
+ tikz("normal-dist.tex", width = 8, height = 4,
+ standAlone = TRUE,
+ packages = c("\\usepackage{tikz}",
+ "\\usepackage[active,tightpage,psfixbb]{preview}",
+ "\\PreviewEnvironment{pgfpicture}",
+ "\\setlength\\PreviewBorder{0pt}",
+ "\\usepackage{amssymb}"))
> u=seq(-3,3,by=.01)
> plot(u,dnorm(u),type="l",axes=FALSE,xlab="",ylab="",col="white")
> axis(1)
> I=which((u<=1))
> polygon(c(u[I],rev(u[I])),c(dnorm(u)[I],rep(0,length(I))),col="red",border=NA)
> lines(u,dnorm(u),lwd=2,col="blue")
> text(-1.5, dnorm(-1.5)+.17, "$\\textcolor{blue}{X\\sim\\mathcal{N}(0,1)}$", cex = 1.5)
> text(1.75, dnorm(1.75)+.25,
+ "$\\textcolor{red}{\\mathbb{P}(X\\leq x)=\\displaystyle{ + \\int_{-\\infty}^x \\varphi(t)dt}}$", cex = 1.5)
> dev.off()

Now we have the graph in another tex file. It is possible to embed the code in a tex file, or to compile the tex file to get a pdf file. I did generate the pdf file.

Here is the tex file I finally get. It is now extremely simple to get your own normal distribution table. Now, I guess it could be possible to use sweave, or knitr. Once I’ll get a copy of Yihui’s book, I’ll try to use it to generate distribution table for my courses !

After four days offline (at least off my blog, see the previous post for more details), I have to face the truth: I am a computer addict. For sure. Here is the diary of the last four days, that were supposed to be without touching my computer, at work and at home. I tried to keep tracks of everytime I had to go on my computer. At home, that was fine (I have decided a few weeks ago that I should not check my email at home, in the evening, and in the morning, so unless I want to read the news, check on Twitter what’s going on, or write a post on my blog, I do not usually spend much time on our computer). But at the office, that’s another story…

• Tuesday, April 2nd

6:17 Wake up, first day of the challenge.

8:17 Time to check if the code used for datascraping (on some websites) did run… It’s not like using my computer. It was just fixing problems for a future research. One minute, just cheking. Well, there was a problem in my code, I have to fix it (it takes much longer than I thought) and then, I ran it again (in order to exact some figures out of almost 200,000 internet pages). One code has been scaping a website for almost a week (and it looks like I have only one third of the data), and the other one, I have to run it everyday, to backup some daily figures from several websites.

8:32 While going to get a coffee, JF shows me the web site of the Antartica Journal of Mathematics (where he was kindly invited to submit an article, and also to apply if he’s willing to join the board) and we go through Rob Hyndman’s warning on his blog, about junk journals. No offense to the webdesigner hired by this journal, but we had a lot of fun on that website… so… nineties.

8:47 I have to go online, on my email account to download a paper I have to review. I had a reminder by an editor this weekend. Print the article. While seeking for the email, I notice that I did receive during the night an invitation to give my opinion about another article, for another journal. Just go briefly through the paper. Even if my (personal) quota for 2013 is already exceeded, I decided to accept to write a referee report. Also check something with a co-author. Quickly. Also read four comments on the blog submited during the night, and approve all of them. Damned! one of them mentioned a preprint related so something discussed in a post. Try to avoid to read the preprint. Go offline. Definitively. Turn on the music. Sacred music, namely In Seculum Longum.

9:39 After alsmost 30 minutes reading articles, I give up. I open LaTeX to type corrections for a couple of chapters I have to write down (and send them as soon as possible, the deadline was last week). So far, no internet. I take a gum, I feel nervous…

9:44 Go online to check how to use enquote{} in latex (it looks like I have some damned babel problems !) Go on several forums. Spent the morning working on my LaTeX files.

12:54 Back from lunch, I turn off the sacred music, swith to Sandinista, by The Clash,

13:14 JF want to check with me the schedule for my graduate courses for Winter 2014 session, one on extreme values and copulas, and one on time series. Have to go on the website of the university, to find what has been planed. Spend also some time seeking for old emails, since the information is only partially online.

14:02 Have to go online, one more time, to buy a ticket for the French railway, for a colleague of mine.

14:17 Mathieu asks me go check the emails received this morning, and to send an email in order to confirm that I will join a meeting about the organisation of a workshop.

14:40 Discussion with Mathieu about resubmiting an article. We need to go on my email account to check several (interesting) comments from the referees.

15:03 This time, I have to launch the internet browser. I have to find out how to avoid fist names in a bibliography, using a bib file, and get only initials. Some more time on LaTeX forums.

15:08 Remember that I have to book a room for a colleague, visiting me this August. Need to check the price on the webiste, and send emails.

15:19 Wrote a recommendation for a former student of mine. Send it.

15:26 Discussion by email with two co-authors, since we did plan to work on a joint paper this month.

15:48 Time to work on slides for a conference at the end of this month, and write codes to produce some graphs. So, still on my computer. But not (really) online.

Finally, I went back home early, to cook, and spend some time with the kids.

• Wednesday, April 3rd

Second day.

8:17 Check again my R code. Run again one of them. Looks like there was a problem.

8:50 Go to give my course. Until noon, I spent the morning producing codes, to show how to compute chain ladder estimates, and explain roots of bootstrap in regression models.

13:03 Buzy red light on my phone when I get back to my office: it looks like some people at the faculty tried to reach me. Usually I do not answer the phone. I hate my phone: you want to reach me, send me an email. OK, in my mailbox, there are a couple of emails from the faculty “please, call us back, regarding the conference you plan to organize”. Damned, how can I tell them I am on a mission, that I try to avoid using my email (and that I have to deal with my phone-phobia at the same time).

13:16 Go online to check code for an R package I want to use to produce nice graphs

13:30 Check quickly my emails, delete 90% of them, answer to one, postpone others. Decide to go to work in a coffee shop, the whole afternoon. While taking the elevator, I started to discuss with a colleague. Looks like I missed an email about a meeting that will take place in the afternoon. I want to work on my chapters. I go to the coffee shop. Nothing serious in the meeting, as far as I understood. I spent the afternoon reading, checking typos in chapters that should be sent soon, and articles on advanced methodes in finance, based on trees… No computer for a few hours ! I did it !

17:49 I have to go online at home, time to pay my HydroQuébec bill. Damned !

17:56 While I am about to log off, I received an email from a former student of mine, with a link to a nice article (entitled “l’informatique, ça s’apprend“) I really want to share it. But I can’t. Less than 36 hours after my mission, that would be a defeat ! Do not go on Twitter !

18:02 Cooking for the kids. Will miss the Montreal Hackathon organized by R users. Wednesdy evening is not a great day to join those social events (could geek meetings be called social ?). My wife has a late course in the evening, and I still try to see how deeply addicted I can be. I finally decided that I will check my emails twice a day. But just to remove spams (or messages like “I am a student in a engineering program from India and I would like to start a PhD with you“, or “the back door in one of the building will be locked during the week-end from 22:15 till 23:42 for security reasons“) and to see if there might be important ones. Probably have to answer some of them, but I’ll try to postpone for most of them.

18:41 started a game of kid audio with the girls. Almost became a pugilat when they started to argue about kettledrum and snare-drum, I wanted to show them on youtube, but finally I gave up (there is an old saying: never interfere in a girl-and-girl fight). Decided to ask the elder to read a story to the youngest, while I was washing the dishes (which is usually the perfect timing for a dvd). Meanwhile, it looks like my son went online, for his music assignment: his teacher is using online videos to help them practicing. Argggg.

• Thursday, April 4th

Third day.

8:29 As usual, checking the code for datascraping… reload the one that did crash (again) during the night

Error in substr(html, 8, 12) :
invalid multibyte input string at '<e9>lair,'

Damned. Moving around 200,000 pages without being caught is difficult. Have to play some music. Gonzales, piano solo.

8:47 Have to check quickly my emails. The problem is that, on average, per week day, I have a bit more than 100 emails (excluding official spam). If I do now scan them, I end up with on thousand emails very quickly… Need to moderate a comment on the blog.

8:50 Still online, checking my emails, bad news about fundings for a student of mine, have to send a couple of emails to find a backup solution.

8:54 Quick discussion by email about copyrights for a chapter in a book

8:54 Have to send also emails to book a room for a colleague who will visit me in August

8:55 Postpone a Skype discussion with a co-author, still trying to avoid unecessary use of my computer.

8:57 Answer an email to schudele a meeting because a student asked to get a grade revision, and an adhoc committee is necessary. Looks like I am part of it.

9:17 Start to write recommendation letters for a Christophe who’s applying for positions in several universities, in France.

9:34 Back on the slides and the R code. On my computer, but offline.

10:11 Email, brief answer to a former student of mine, who might be interesting to share some datasets, but it looks like there might some confidentiality issues. I wanted to work on those data with a student, in September. Have to find an alternative.

10:45 Discussion with student in master program (face to face this time). Have to go on Dropbox to download a pdf file he wrote, and to download a couple of paper to check the proofs.

13:05 Work with Amadou, my phd student, need to find a pdf version of a book, since the property is clearly given there, but the book is out of print (no way to get a paper copy). Also go online to find a reference on a complicated model.

14:40 Email from Fred, about a reserving technique that seems new, from a paper he just discovered.

15:03 Upload on slideshare some (old) slides that do the same thing that this new paper (aren’t there anyone checking before publishing papers that results are really new ?)

15:14 Play Rodrigo y Gabriela, need something punchy to finish my day

16:56 Received an email from the immigration department, and I have to go to their website to find a doctor for a medical examination of the whole family.

19:26 Request by email from financial services to get the exact amount (in Euros) taken from my credit card. Have to go on my bank account, online.

19:43 Disccussion by email with Frédéric about one year uncertainty and bootstrap with overdispersed Poisson models

• Friday, April 5th

Last day of the test. Fourth day.

06:18 My daughter wakes up and tells me it is unfair to have snow for her birthday. Have to go on http://meteomedia.com/weather/… to check for weather forecasts. Hopefully, we should have nice weather this afternoon…

08:16 Once again, checking R codes, still running this time ! Great ! Time to play some music.  Air, Premiers Symptômes.

8:38 Email regarding next week jury, checking legal aspects

9:04 I have to print bank informations, that I did download yesterday, check orders placed on Amazon in the last 3 months, scan documents, send them to financial services.

9:13 Have to check for an account number with financial services

9:15 Brief email to some contributors of a book that I should edit this year.

9:21 Lauch Skype, have to discuss with coauthors in France.

13:10 Update the syllabus for the course I will give in September. Decided to write that cellphone will not be allowed during the class.

14:12 Work with Ben, a master student, on a paper. Need to scan notes I have written down to send him by email (this time, we did work together without using my computer).

14:38 Finished my recommendation letters for Christophe (for positions in France, more than a dozen recommendations). Have to send them individualy by email. Have also to find some email adresses that I do not have.

14:49 Received an email claiming that I cannot give my (graduate) course on extreme value in 2014, not in the official programm. Have to spend some time checking why. It seems that it has been removed from the list, and that the code has changed (hopefully JF was online to check that information much more efficiently than I would).

14:59 Received an answer from one of the colleagues I just sent recommendation to. We used to be students at the same time, a few years ago. Write a short email back to give (personal) news.

16:19 Started to type sketches of what can be the final exam for my course on GLM for actuarial science. So far 5 questions. Need to find about 40…. Will take a while.

Mission aborted.

I finally left the office later on, to pick up my son and bring him to his fencing course (to prepare for the Jeux de Montréal that will take place tomorrow). Went back home, then, for a cake with candles for my daughter’s birthday. Later on, had to spend some time online, on the blog, for my students. And started to type this post. So, here was the story of the past four days. I have to admit that looking back at those four days is quite informative :

1. I cannot work without a computer, I can hardly work offline. No only for my research, to get help from forums on R and LaTeX, or to seek articles (the time I spent in Paris at the library making photocopies of old articles is clearly over).
2. I do not only need a computer to write R codes, and produce LaTeX slides and artices. I need a computer… for everything. To plan meetings, for social interactions with colleagues, to find some help, to find a theorem, to book an hotel, etc.
3. I understand more clearly why I am so unproductive in terms of research ! it looks like I spend (I was about to say waste) a lot of time on administrative tasks. Small tasks. But adding a lot of small tasks, finally, it is difficult to have 4 or 5 consecutive hours to work exclusively on some research project.