Tag Archives: computer

The law of small numbers

In insurance, the law of large numbers (named loi des grands nombres initially by Siméon Poisson, see e.g. http://en.wikipedia.org/…) is usually mentioned to legitimate large portfolios, because of pooling and diversification: the larger the pool, the more ‘predictable’ the losses will be (in a given period). Of course, under standard statistical assumption, namely finite expected value, and independence (see http://freakonometrics.blog.free.fr/…. for a discussion, in French). Since in insurance, catastrophes are usually rare – and extremely costly – and actuaries might be interested to model occurrence of that small number of events (see e.g. Aldous’ book on that specific topic, that can be downloaded from http://stat.berkeley.edu/…). The theorem behind is sometimes called the law of small numbers (from the book published by Ladislaus Bortkiewicz, but we’ll get back to that story later on, see also Whitaker (1914) http://biomet.oxfordjournals.org/… or the book recently published by Michael Falk, Jürg Hüsler and Rolf-Dieter Reiss).

  • The Poisson distribution

The so-called Poisson distribution (see http://en.wikipedia.org/…) was introduced by Siméon Poisson in 1837 (in Recherches sur la Probabilité des Jugements en Matière Criminelle et en Matière Civile, Précédées des Règles Générales du Calcul des Probabilités, see http://gallica.bnf.fr/…). But it had been defined more than a century before, by Abraham De Moivre, in 1711, in De Mensura Sortis seu; de Probabilitate Eventuum in Ludis a Casu Fortuito Pendentibus (see e.g. the review in http://www.jstor.org/…). Let https://latex.codecogs.com/gif.latex?N denote a counting random variable, then it said to be Poisson distributed if there is https://latex.codecogs.com/gif.latex?\lambda\in(0,\infty) such that


De Moivre obtained that distribution from an approximation of the binomial distribution. Recall that the binomial distribution is a standard distribution in actuarial science, for instance to model the number of deaths among https://latex.codecogs.com/gif.latex?n insured. If individual death probabilities are identical, say https://latex.codecogs.com/gif.latex?p, and if deaths are independent events, then

And if https://latex.codecogs.com/gif.latex?n\rightarrow\infty and  https://latex.codecogs.com/gif.latex?np\rightarrow%20\lambda, then

https://latex.codecogs.com/gif.latex?\mathbb{P}(N=k)\rightarrow%20e^{-\lambda}\frac{\lambda^k}{k!}Again, this is an asymptotic theorem, which is valid when we have a lot of observations (https://latex.codecogs.com/gif.latex?n\rightarrow\infty), but also, the probability of occurrence should be extremely small (since https://latex.codecogs.com/gif.latex?p\sim\lambda/n), which is why to use the term small numbers. Siméon Poisson was not interested by mathematical approximations: his main point was to get a distribution with nice goodness of fit properties for the data he was working on. He wanted to get a better understanding of cours d’assises (jury panel, might be a valid translation of the French term). A jury consists of 12 jurors who voted to determine whether a defendant was guilty. When guilt was predominant, with at least 8 votes against 4, the defendant was convicted (which was 47% of criminal cases). 5 with 7 votes against, the opinion of professional judges was requested (11% of criminal trials again). Using these statistics we can demonstrate that a defendant brought before an assize court is guilty of the order of 68%, and the probability that a juror is not wrong by voting (condemning an innocent or releasing a culprit) was about 54%. He sought to calculate the probability that a defendant is wrongfully convicted, and gets 2%. And 28% of exonerated defendants are in fact guilty. Siméon Poisson introduced this law to get probabilities easily. But the law he considered is central in probability….

  • The law of small numbers

The heuristic of the main theorem, related to the Poisson distribution is the following: let https://latex.codecogs.com/gif.latex?X_1,%20\cdots,X_n denote i.i.d random variables taking values in https://latex.codecogs.com/gif.latex?%20\mathbb{R}^d (in a general setting, one component can be the time, the other one an upper region of interest, where some stochastic process might be). Let https://latex.codecogs.com/gif.latex?\mathcal{A}_n\subset\mathbb{R}^d. If  https://latex.codecogs.com/gif.latex?\mathbb{P}(X_i%20\in%20\mathcal{A}_n)\rightarrow%200 as https://latex.codecogs.com/gif.latex?n\rightarrow\infty (or https://latex.codecogs.com/gif.latex?\mathbb{P}(X_i%20\in%20\mathcal{A}_n)=O(n^{-1}) to be a little bit more specific about the assumptions), let https://latex.codecogs.com/gif.latex?N denote the (random variable characterizing) count of events https://latex.codecogs.com/gif.latex?\{X_i%20\in%20\mathcal{A}_n\}, then https://latex.codecogs.com/gif.latex?N can be approximated by a Poisson distribution with parameter https://latex.codecogs.com/gif.latex?\lambda%20=n%20\times%20\mathbb%20P(X_i%20\in%20\mathcal{A}_n).
The heuristic is that if we consider a large number of observations, and if we count how many are in a given (small) region, then the number of such observations is Poisson distributed.


If we run some simulations,

>  n=1000
>  ns=100000
>  N=rep(NA,ns)
> for(s in 1:ns){
+ X=runif(n)*10-1.5
+ Y=runif(n)*10-1.5
+ I=(X^2+Y^2)<1
+ N[s]=sum(I)
+ }
> hist(N,breaks=0:60,probability=TRUE,col="yellow")
> mean(N)
[1] 31.41257

The parameter of the Poisson distribution is the area of the yellow disk, over the area of the square, i.e.

> (lambda=10*pi)
[1] 31.41593
> lines(0:60-.5,dpois(0:60,lambda),type="b",col="red")


To get an interpretation related to insurance modeling, let https://latex.codecogs.com/gif.latex?\mathcal{A} denote an upper layer in a reinsurance contract, i.e. https://latex.codecogs.com/gif.latex?\mathcal{A}=\{x%3Ed\} for some deductible https://latex.codecogs.com/gif.latex?d. Let https://latex.codecogs.com/gif.latex?X_i‘s denote individual losses. Then the number of claims that hit this upper layer can be modeled with a Poisson distribution. More precisely, if deductible https://latex.codecogs.com/gif.latex?d becomes extremely large (and https://latex.codecogs.com/gif.latex?\mathbb{P}(X_i%20\in%20\mathcal{A})\rightarrow%200), we obtain the point-over-threshold model in extreme value theory (see e.g. http://brale.math.hr/~iugrina/… or  http://fire.nist.gov/bfrlpubs/…): if https://latex.codecogs.com/gif.latex?N has a Poisson distribution and, conditionally on https://latex.codecogs.com/gif.latex?Nhttps://latex.codecogs.com/gif.latex?X_1,\cdots,X_N are independent identically distributed generalized Pareto random variables, then https://latex.codecogs.com/gif.latex?\max\{X_1,\cdots,X_N\} has the generalized extreme value distribution. Thus, exceedances models (for rare events) are closely related to Poisson processes.

  • The Poisson process

As mentioned above, the Poisson distribution appears when events occur somehow randomly and independently, over time. It is then natural to study the time between two occurences (or two claims, in an insurance context).

  • Poisson distribution, and claims occurrence

It is neither Siméon Poisson nor De Moivre, but Ladislaus Von Bortkiewicz who first mentioned the Poisson distribution as the law of small numbers. In 1898 (see https://archive.org/…), he studied the number number of soldiers killed by being kicked by a horse, from 1875 till 1894, in 200 corps (more precisely 10 corps over 20 ans).

He did obtain the following distribution (here, the parameter of the Poisson distribution is 0.61, i.e. the average number of death per year)

number of
death per
0 109 108.67
1 65 66.21
2 22 20.22
3 3 4.11
4 1 0.63
5 and more 0 0.08

It is possible to find a lot of cases where the Poisson distribution fits extremely well. For instance, if we consider the number of hurricanes, that landed in Florida after 1850,

number of
per year
0 30 27.16
1 48 47.99
2 37 42.41
3 29 24.98
4 8 11.03
5 3 3.90
6 3 1.15
7 1 0.29
8 and more 0 0.08
  • Poisson distribution, and return period

The return period was introduced by Emil Gumbel, in hydrology, to link probabilities and durations (see e.g. http://freakonometrics.blog.free.fr/…). A decennial event has an occurence probability of 1/10. 10 is then the average waiting time before occurence. This does not mean that the event will not occur before 10 years, or has to occur before 10 years. Consider a return period https://latex.codecogs.com/gif.latex?T (in years), then the yearly probability of non-occurrence is https://latex.codecogs.com/gif.latex?1-(1/T).

And the probability of non-occurence over https://latex.codecogs.com/gif.latex?n years is then https://latex.codecogs.com/gif.latex?1-[1-(1/T)]^n. It is standard to summarize this property with the following table,

return period https://latex.codecogs.com/gif.latex?T

Number of years (https://latex.codecogs.com/gif.latex?n) without catastrophes

10 20 50 100 200
10 65.1% 40.1% 18.3% 9.6% 4.9%
20 87.8% 64.2% 33.2% 18.2% 9.5%
50 99.5% 92.3% 63.6% 39.5% 22.5%
100 99.9% 99.4% 86.7% 63.4% 39.5%
200 99.9% 99.9% 98.2% 86.6% 63.3%

The diagonal in the table above is extremely interesting. It looks like there is some kind of convergence towards a limiting value (here  63.2%). Indeed, the number of events observed over n years have a Binomial distribution, with probability https://latex.codecogs.com/gif.latex?1/T=1/n, which will converge towards the Poisson distribution with parameter 1. The probability of not having a catastrophe is then https://latex.codecogs.com/gif.latex?1-\exp(-1), which is equal to 0.632.

  • Rare probabilities and the Poisson distribution

The Poisson distribution keeps appearing when computing probabilies of rare events. For instance, the probability to have at least one incident in a nuclear plant in France, over a 50 year period. Assume that the annual probability of an incident in a reactor https://latex.codecogs.com/gif.latex?p is small, e.g. 0.05%. Assume further that reactors are independent among them, and in time. The probability to have an incident over 80 reactors in 50 years is (exactly)


Of course, a linear approximation is not correct (even if it was mentioned in some French newspaper, as explained in an old post http://freakonometrics.blog.free.fr/…)


On the other hand

https://latex.codecogs.com/gif.latex?\mathbb%20P(N\neq 0)=1-(1-p)^{50\times80%20}%20\sim1-\exp\left(-50\times80\times%20p%20\right)

> p=0.00005
> 1-(1-p)^(50*80)
[1] 0.1812733
> 1-exp(-50*80*p)
[1] 0.1812692

which is the probability that https://latex.codecogs.com/gif.latex?N is null when https://latex.codecogs.com/gif.latex?N has a Poisson distribution with parameter https://latex.codecogs.com/gif.latex?\lambda=50\times80\times%20p. We clearly see here an application of De Moivre’s approximation in risk management.

Another way of looking at this problem is based on the following idea: given the fact that in 45 years of observations on 450 reactors worldwide (roughly), three major accidents were observed including Three Mile Island (1979) and Fukushima (2011), i.e. the average time between accidents can be estimated at 16 years. For a single reactor, we can assume that the average time to wait before an incident is 450 times 16 years, i.e 7200 years. Or the probability to have one incident, over one year, for one reactor is 1 over 7200 (this is the idea behind the return period concept). If we assume that the arrival of accidents occurs randomly and independently of each other (as defined above) then the number of major accidents observed over a period of 50 years in France follows a Poisson distribution with parameter 50 / (7200/80). Also, the probability of having no major accident over 50 years, with 80 reactors can be estimated by



> 1-exp(-50*80/7200)
[1] 0.4262466

(keeping in mind all the uncertainty around the estimated waiting time before a major accident to a single reactor!).


I will be in Amsterdam for the end of this week. I will be in the jury of the PhD defense of Julien Tomas, entitled “Quantifying Biometric Life Insurance Risks With Non-Parametric Smoothing Methods” (the thesis will probably be online soon). But before, I will give a talk at the actuarial seminar at UvA. My visit last time was a real pleasure, and it should be the same this time too. I will give a talk this Thursday on “R for actuarial science“. The slides can be downloaded from here.

De la non-connexité du Vaucluse

Avant-hier, José, un ancien collègue rennais me faisait noter une bizarrerie de la cartographie (et me posait des questions sur l’impact sur les cartes faites avec R). En fait, il m’a fait découvrir que le département du Vaucluse n’était pas connexe. Comme on le voit sur la carte ci-contre, il y a l’enclave des papes, qui est enclavé dans la Drome, mais administrativement rattachée au Vaucluse. Étonnant non ?

Maintenant avec R, ce genre de choses existent. Par exemple, il est possible de travailler avec les îles, qui sont rattachées à tel ou tel département. Regardons ce qui se passe ici, avec les cartes standards de R,

>  library(maps)
>  france = map(database="france")
>  france$names
  [1] "Nord"                                
  [2] "Pas-de-Calais"


 [92] "Gard"                                
 [93] "Vaucluse"                            
 [94] "Tarn-et-Garonne"                     
 [95] "Alpes-Maritimes"                     
 [96] "Vaucluse"                            
 [97] "Tarn"                                
[108] "Hautes-Pyrenees"                     
[109] "Var:Iles d'Hyeres:I. du Levant"      
[110] "Var:Iles d'Hyeres:I. de Porquerolles"
[111] "Var:Iles d'Hyeres:I. de Port Cros"   
[112] "Haute-Corse"                         
[113] "Pyrenees-Orientales"                 
[114] "Corse du Sud"

On voit que le Vaucluse apparaît deux fois dans la liste des départements. Pour les îles, elles sont rattachées à un département avec un nom spécifique (comme on le voit sur l’île de Porquerolles, par exemple). Mais pas l’enclave des papes. En fait, si on cherche le Vaucluse, il apparaît deux fois

>  which(substr(tolower(france$names),1,5)=="vaucl")
[1] 93 96

Aussi, si on colore le Vaucluse, c’est le département tout entier (avec l’enclave) qui ressort,

 Le code est ici

>  dpt="Vaucluse"
>  couleur="red"
>  match=match.map(france,dpt)
>  color=couleur[match]
>  map(database="france", fill=TRUE, col=color)

On peut aussi faire ressortir l’enclave. Pour cela, il suffit d’aller demander de colorer de manières différentes les deux régions,

>  match[which(match==1)[2]]=2
>  couleur=c("blue","red")
>  color=couleur[match]
>  map(database="france", fill=TRUE, col=color)

Ah, la joie des cartes avec R…


Simulation de séries temporelles

Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons  un processus autorégressif d’ordre 1,  où  est un bruit blanc, stationnaire, i.e.  appartient à l’intervalle . Le code pour simuler un tel processus est

phi1= .85
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}

ou avec un autocorrélation au premier ordre négative,

phi1= -0.7

On peut aussi regarder un processus autorégressif au second ordre,

sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres).

phi1=  0.3
phi2=  0.5
for(t in 3:n){

Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre,  où  est un paramètre dans .

theta1=  .8
for(t in 2:n){

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
for(t in 3:n){


Inference and autoregressive processes

Consider a (stationary) autoregressive process, say of order 2,

for some white noise  with variance . Here is a code to generate such a process,

> phi1=.5
> phi2=-.4
> sigma=1.5
> set.seed(1)
> n=240
> WN=rnorm(n,sd=sigma)
> Z=rep(NA,n)
> Z[1:2]=rnorm(2,0,1)
> for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}

Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.

  • using least square regression

A natural idea is to see here a regression model, and thus, 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)

lm(formula = Y ~ 0 + X1 + X2, data = base)

Min      1Q  Median      3Q     Max
-4.3491 -0.8890 -0.0762  0.9601  3.6105

Estimate Std. Error t value Pr(>|t|)
X1  0.45107    0.05924   7.615 6.34e-13 ***
X2 -0.41454    0.05924  -6.998 2.67e-11 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.449 on 236 degrees of freedom
Multiple R-squared: 0.2561,	Adjusted R-squared: 0.2497
F-statistic: 40.61 on 2 and 236 DF,  p-value: 6.949e-16

> regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma
[1] 1.449276
  • 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

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,]  0.4517579
[2,] -0.4155920

Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)

> estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.445706

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.

  • using (conditional) likelihood estimators

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

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)
[1]  0.4509685 -0.4144938  1.4430930

[1] 425.0164

function gradient
88       NA

[1] 0


Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient

for the second one,

and finally for the standard deviation of the innovation process

All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).

Les députés sont-ils à l’image de la population

Beaucoup de choses ont été écrites sur le fait que les députés ne sont pas vraiment le reflet de la population, que ce soit en terme de profession, de sexe, d’origine, d’age, etc. La liste pourrait être longue. Il y a plusieurs mois, j’avais commencé à regarder le profil des députés, par age. En effet, le site http://assemblee-nationale.fr/ permet d’accéder à des données sur tous les députés, depuis la Révolution. Y compris leur date de naissance. En croisant ces données avec des données de population, par exemple via http://www.mortality.org/, on peut comparer la répartition des ages des députés, avec la répartition des ages de la population.

Pour les amateurs, le code pour récupérer les données (ou au moins les dates de naissance des députés) ressemble à



Maintenant, pour être tout à fait honnête, je ne suis pas certain de ce qui est vraiment renvoyé, et j’ai des doutes que cela correspondent réellement à la requête faite. En effet, même si je demande à avoir la liste des députés après l’élection, j’ai trop de monde… mais peut-être est-ce du aux décès éventuels, et il est possible que l’ensemble des députés qui ont siégé pendant la mandature apparaissent dans le résultat de la requête.

Sur la figure suivante on voit, sur plusieurs élections depuis plus de 100 ans, comment les deux distributions se déforment, avec en rouge la distribution de l’age des députés, et en bleu, la distribution de la population française, dans son ensemble (population de plus de 18 ans)

Si on veut tout suivre sur un graphique, au lieu de se regarder une animation, on peut représenter les différents quantiles (10%, 25%, 75% et 90%, retenus sur la population de plus de 18 ans, et l’age médian, au centre), avec la population française l’année de l’élection,

et l’ensemble des élus au parlement,

Si on veut faciliter la comparaison, on peut se contenter de visualiser l’évolution des ages moyens,

ou encore, du ratio (en % de différence) entre l’age moyen des députés, et celui de l’ensemble de la population.

Sur ce graphique, on voit que depuis 30 ans, l’age moyen des députés croit plus vite que celui de la population: la population français vieilli, mais moins que ses députés… La gérontocratie perdure donc en France. En espérant que cela ne débouche pas sur le clash générationnel que l’on semble observer ces temps-ci au Québec…

Finding Waldo, a flag on the moon and multiple choice tests, with R

I have to admit, first, that finding Waldo has been a difficult task. And I did not succeed. Neither could I correctly spot his shirt (because actually, it was what I was looking for). You know, that red-and-white striped shirt. I guess it should have been possible to look for Waldo’s face (assuming that his face does not change) but I still have problems with size factor (and resolution issues too). The problem is not that simple. At thehttp://mlsp2009.conwiz.dk/ conference, a price was offered for writing an algorithm in Matlab. And one can even find Mathematica codes online. But most of the those algorithms are based on the idea that we look for similarities with Waldo’s face, as described in problem 3 on http://www1.cs.columbia.edu/~blake/‘s webpage. You can find papers on that problem, e.g. Friencly & Kwan (2009) (based on statistical techniques, but Waldo is here a pretext to discuss other issues actually), or more recently (but more complex) Garg et al. (2011) on matching people in images of crowds.

What about codes in R ? On http://stackoverflow.com/, some ideas can be found (and thank Robert Hijmans for his help on his package). So let us try here to do something, on our own. Consider the following picture,

With the following code (based on the following file) it is possible to import the picture, and to extract the colors (based on an RGB decomposition),

> library(raster)
> waldo=brick(system.file("DepartmentStoreW.grd",
+ package="raster"))
> waldo
class       : RasterBrick
dimensions  : 768, 1024, 786432, 3 (nrow,ncol,ncell,nlayer)
resolution  : 1, 1  (x, y)
extent      : 0, 1024, 0, 768  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values      : C:\R\win-library\raster\DepartmentStoreW.grd
min values  : 0 0 0
max values  : 255 255 255

My strategy is simple: try to spot areas with white and red stripes (horizontal stripes). Note that here, I ran the code on a Windows machine, the package is not working well on Mac. In order to get a better understanding of what could be done, let us start with something much more simple. Like the picture below, with Waldo (and Waldo only). Here, it is possible to extract the three colors (red, green and blue),

> plot(waldo,useRaster=FALSE)

It is possible to extract the red zones (already on the graph above, since red is a primary color), as well as the white ones (green zones on the graphs means a white region on the picture, on the left)

# white component
white = min(waldo[[1]] , waldo[[2]] , waldo[[3]])>220
focalswhite = focal(white, w=3, fun=mean)

# red component
red = (waldo[[1]]>150)&(max(  waldo[[2]] , waldo[[3]])<90)
focalsred = focal(red, w=3, fun=mean)

i.e. here we have the graphs below, with the white regions, and the red ones,

From those two parts, it has been possible to extract the red-and-white stripes from the picture, i.e. some regions that were red above, and white below (or the reverse),

# striped component
striped = red; n=length(values(striped)); h=5
focalsstriped = focal(striped, w=3, fun=mean)

So here, we can easily spot Waldo, i.e. the guy with the red-white stripes (with two different sets of thresholds for the RGB decomposition)

Let us try somthing slightly more complicated, with a zoom on the large picture of the department store (since, to be honest, I know where Waldo is…).

Here again, we can spot the white part (on the left) and the red one (on the right), with some thresholds for the RGB decomposition

Note that we can try to be (much) more selective, playing with threshold. Here, it is not very convincing: I cannot clearly identify the region where Waldo might be (the two graphs below were obtained playing with thresholds)

And if we look at the overall pictures, it is worst. Here are the white zones, and the red ones,

and again, playing with RGB thresholds, I cannot spot Waldo,

Maybe I was a bit optimistic, or ambitious. Let us try something more simple, like finding a flag on the moon. Consider the picture below on the left, and let us see if we can spot an American flag,

Again, on the left, let us identify white areas, and on the right, red ones

Then as before, let us look for horizontal stripes

Waouh, I did it ! That’s one small step for man, one giant leap for R-coders ! Or least for me… So, why might it be interesting to identify areas on pictures ? I mean, I am not Chloe O’Brian, I don’t have to spot flags in a crowd, neither Waldo, nor some terrorists (that might wear striped shirts). This might be fun if you want to give grades for your exams automatically. Consider the two following scans, the template, and a filled copy,

A first step is to identify regions where we expect to find some “red” part (I assume here that students have to use a red pencil). Let us start to check on the template and the filled form if we can identify red areas,

exam = stack("C:\\Users\\exam-blank.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)

First, we have to identify areas where students have to fill the blanks. So in the template, identify black boxes, and get the coordinates (here manually)

exam = stack("C:\\Users\\exam-blank.png")
black = max(  exam[[1]] ,exam[[2]] , exam[[3]])<50
focalsblack = focal(black, w=3, fun=mean)

The blue points above are where we look for students’ answers. Then, we have to define the vector of correct answers,

points(CORRECTX, CORRECTY,pch=16,col="red",cex=1.3)

Now, let us get back on red areas in the form filled by the student, identified earlier,

exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=5, fun=mean)

Here, we simply have to compare what the student answered with areas where we expect to find some red in,

points(CORRECTX, CORRECTY,pch=1,

Crosses on the graph on the right below are the answers identified as correct (here 13),

> icorrect=values(red)[(750-CORRECTY)*
+ 610+(CORRECTX)]
> points(CORRECTX[icorrect], CORRECTY[icorrect],
+ pch=4,col="black",cex=1.5,lwd=1.5)
> sum(icorrect)
[1] 13

In the case there are negative points for non-correct answer, we can count how many incorrect answers we had. Here 4.

> iuncorrect=values(red)[(750-UNCORRECTY)*610+
> sum(iuncorrect)
[1] 4

So I have not been able to find Waldo, but I least, that will probably save me hours next time I have to mark exams…

Correlations, dimension, and risk measure

Yesterday, while I was attending the IFM2 conference, at HEC Montreal, I heard a nice talk about credit risk, and a comparison between contagion (or at least default correlation), for corporate and retail companies (in the US). And it was mentioned that default correlation was much lower for retail companies than it could be for corporate risk. In a discussion that followed those slides, it was mentioned that banks in the US should actually have been working more with those small firms, since contagion risk was much lower.

A problem here is that the link between correlation, risk and dimension is rather complicated:

  • corporate means a small number of firms, high correlation (and possible large individual losses)
  • retail means a large number of firms (even perhaps extremely large), lower correlation (and small individual losses)

A simple model for default models is based on the assumption that we deal with an exchangeable portfolio (as in a previous post). With the following code, given an (individual) default probability, a default correlation, and a number of firms, it is possible to calculate the probability to have more than a given number of defaults.

 stop.on.error =  FALSE)$value}

a=m*(1-r)/r ;
 for(i in 0:n){
 return(sum(V[1:(x+1)])) }

It is possible to calculate, for a large range of correlations, the probability to have – at least – 20% of default in the portfolio (in order to compare things that are comparable).

for(i in 1:length(A)){

With 20 firms (corporate) we want to have at least 4 defaults, while with 1000 firms (retail) there should be 200 defaults. As mentioned in the previous post, the relationship between correlation and quantiles of sums is not simple. Hence, it might not be monotone. The dotted line is the probability to have at least 4 defaults when default correlation is 50% (around 10%). The plain line is the probability to have at least 200 defaults, as a function of the correlation,


In that case, with only a correlation of 10% among retail firms, the probability of having 20% defaults is the same as the same probability for corporate, but with 50% correlation… One should remember that in portfolio analysis, the links between correlation, dimension and risk measure is a sensitive issue…

Régression médiane et géométrie

Suite à mon exposé d’hier, Pierre me demandait d’argumenter un point que j’avais évoqué oralement sans le justifier: “ une régression médiane passe forcément par deux points du nuage “. Par exemple,


Essayons de le justifier… de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe “entre” les points (on admettra que la régression médiane sépare l’espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).


c’est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c’est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

for(i in seq(-2,3,by=.5)){

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))

Formellement, l’optimum est ici

> optimize(d,lower=-5,upper=5)
[1] -0.2000091

[1] 2.200009

Bref, retenons cette courbe, et notons qu’elle passe par un des points,

Et c’est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l’on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu’il n’y a qu’un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s’arrêter dès qu’on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

for(i in seq(-1,5,by=.25)){

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
[1] 2.500018

[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,


Pour comprendre pourquoi, comparons les deux cas,

  • en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue

  • en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c’est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu’à l’atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c’est toujours pareil…

Simulation d’un processus de Lévy, et discrétisation

Avec @renaudjf, on discutait l’autre jour de la simulation d’un processus de Lévy. Et on se posait la question d’un algorithme optimal pour combiner un processus de Poisson (ou un process Poisson composé) avec un processus de Wiener (avec éventuellement un drift, voire une diffusion plus générale). En fait, pour générer des processus de Poisson, j’ai toujours eu l’habitude de simuler les durées entre sauts (avec des lois exponentielles, indépendantes, comme dans des vieux billets). Jean François me suggérait d’utiliser une propriété d’uniformité des sauts sur un intervalle de temps donné, conditionnellement aux nombres de sauts.

Commençons par la première piste. On peut générer un processus de Wiener, éventuellement avec un drift, et à coté, on peut générer les lois exponentielles  (qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici

où . On commence par générer  en notant que

où les incréments  sont Gaussiens (centrés et de variance ) et indépendants les uns des autres. Quant aux durées entre sauts, les  , ce sont des lois exponentielles indépendantes, de moyenne . Voilà le code qui permet de générer les trois composantes,


Le hic est que pour le processus de Wiener, on a du discrétiser, alors que pour le processus de Poisson composé, non. Pourtant, il va bien se ramener sur une même échelle de temps. Une première piste est de créer vraiment la fonction 


et pour faire un dessin ensuite, c’est un jeu d’enfant


Une autre possibilité est d’utiliser une propriété d’uniformité du processus de Poisson que j’évoquais en introduction. Car le processus de Poisson vérifie une propriété remarquable: si  est la date où survient le ième saut, , alors conditionnellement au fait que , les variables  correspondent aux statistiques d’ordres de  variables indépendantes, uniformément distribuées sur , i.e.

Cette propriété se trouve dans Wolff (1982). L’idée de la démonstration est relativement simple. Commençons avec un (unique) saut, alors pour ,

i.e. on retrouve la fonction de répartition d’une loi uniforme sur . On itère ensuite avec 2 sauts, 3 sauts, etc.

La traduction en R de cette idée est tout simplement (car on se place sur )


Ensuite, une stratégie est de discrétiser le processus de Poisson, avec le même pas de temps que le processus de Wiener,


On retrouve la même trajectoire qu’auparavant



Sauf qu’on a eu de la chance. Avec cette procédure, il ne faut pas que l’on ait deux sauts dans le même intervalle de temps ! Bon, il est vrai qu’une caractérisation du processus de Poisson est que

et donc on doit avoir très peu de chance d’avoir deux sauts au même instant d’autant plus que le pas de temps est petit. Mais “peu de chance” ne veut pas dire nul, et si on génère des milliers de trajectoires, la probabilité d’avoir une fois un soucis n’est pas négligeable.

Jean-François a eu l’idée brillante de proposer de tirer non pas des lois uniformes sur , mais des lois uniformes discrètes, dans


sans remise afin d’éviter d’avoir deux sauts au même moment. L’idée était séduisante, mais il me restait à etre convaincu que les durées entre sauts… continuaient à être (presque) exponentielles.

Pour ça, on peut faire quelques tests. Par exemple, générer quelques simulations pour avoir une centaines de sauts (et donc une centaine de durées entre sauts), puis faire un test de loi exponentielle (de moyenne )

for(ns in 1:20){

On fait ici 20 boucles car on avait fixé


et j’avais dit que je voulais une centaine d’observations pour faire un test de loi (ce qui est purement arbitraire, on en conviendra). On peut ensuite faire un test d’ajustement de loi exponentielle,


Si on répète un grand nombre de fois, en changeant le pas de temps (ou le nombre de subdivisons de l’intervalle de temps), on notera qu’effectivement, avec un grand pas de temps (à gauche ci-dessous) on rejette souvent, voire presque toujours, l’hypothèse de loi exponentielle. Mais qu’assez vite, c’est une hypothèse qui n’est pas invraisemblable,

Je ne sais pas entre ces deux fonctions ce qui serait le plus rapide, mais on a deux jolis algorithmes pour générer les processus de Lévy. Non ?

Sondages et prévisions de séries temporelles

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant


(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
for(k in 2:4){

Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste


Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

for(k in 2:4){

avec, là encore, la commande suivante pour lisser


Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),


On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

for(k in c(1,3:4)){

Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via


On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)


On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

variance= -proba%*%t(proba)
for(i in 1:nrow(M)){
for(k in 1:4){

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

for(k in 2:4){

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

for(k in 2:4){

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes


Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

ETS = ets(X)
for(k in 2:4){

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

for(sim in 1:20){

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

for(sim in 1:1000){

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J’avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu’on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n’est pas forcément l’unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c’est d’ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), “when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.”

N’étant toujours pas convaincu, j’évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.


https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")


Sur la figure ci-dessus, non visualise l’impact de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex14.gif, pour un horizon https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse12.gif est alors https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1


> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c’est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.

HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])

Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE

a 805.0389

> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.




Dans ce cas, la prédiction est https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse13.gif. Le code pour faire un lissage double est là encore assez simple,

> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }

Sur la figure suivante, on visualise l’évolution du lissage en fonction de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse15.gif (https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif étant ici fixé),
(le lissage simple – avec le même poids https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Visualisation des résidus pour des données spatiales

Mardi, nous allons travailler un peu sur la modélisation du nombre d’homicides aux États-Unis, à partir de la base

> US=read.table("http://freakonometrics.free.fr/US.txt", 
+ header=TRUE,sep=";")

(je renvoie sur le précédant billet pour un descriptif précis). Idéalement, ça serait parfait si on pouvait visualiser sur une carte les variables. Pour cela, il faut rajouter une colonne à notre base, avec le nom complet des états,

> abreviation=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/etatus.csv", 
> header=TRUE,sep=",") 
> US$USPS=rownames(US) 
> US=merge(US,abreviation) 
> US$nom=tolower(US$NOM)

Cette fois, on va pouvoir faire de la cartographie, les noms de nos états étant (presque) les mêmes que ceux des cartes de R,

> library(maps) 
> VL0=strsplit(map("state")$names,":") 
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,US$nom)

Cette fois-ci, on a toutes les informations pour faire une carte, avec une couleur fonction de la variable d’intérêt (espérance de vie à la naissance, taux d’homicides, taux illettrisme, etc).

> library(RColorBrewer) 
> carte=function(V=US$Murder,titre= 
+ "Taux d'homicides aux Etats-Unis"){ 
+ variable=as.numeric(as.character(cut(V, 
+ quantile(V,seq(0,1,by=1/6)),labels=1:6))) 
+ niveau=variable[ETAT] 
+ couleur=rev(brewer.pal(6,"RdBu")) 
+ noml=levels(cut(V,quantile(V,seq(0,1,by=1/6)))) 
+ map("state", fill = TRUE, col=couleur[niveau]); 
+ legend(-78,34,legend=noml,fill=couleur, + cex=1,bty="n"); 
+ title(titre)}

Commençons par l’analyse du nombre de jours durant lesquels la température passe en dessous de 0°C, par an, en moyenne, dans la capitale (ou la plus grande ville) de l’état, afin de tester notre fonction,

> carte(US$Frost, + titre="Nombre de jours de gel par an")

Pour le taux d’homicide (qui est la variable par défaut) on a

> carte()

Ça sera notre variable d’intérêt lors de la modélisation de mardi. Enfin, on peut lancer une régression, et représenter spatialement les résidus,

> reg=lm(Murder~.-NOM-USPS-nom,data=US) 
> regs=step(reg) 
> carte(residuals(regs), 
+ titre="Résidus de la régression")

Nous voila équipés pour commencer l’économétrie spatiale…

Pour aller maintenant un peu plus loin dans la modélisation, je vais rajouter une variable qualitative, par exemple l’appartenance politique du gouverneur en 1977 (les données datent de cette époque). Les données ont été extraites de wikipedia, suite aux élections de 1974, 1975, 1976 et 1977,

> GV=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/governor.csv", 
+ header=TRUE,sep=";") 
> etat=strsplit(as.character(GV$State),"-") 
> listeetat=rep(NA,nrow(GV)) 
> for(i in 1:nrow(GV)){ 
+ listeetat[i]=etat[[i]][1] 
+ } 
> indice=which(is.na(listeetat)==FALSE) 
> basegv=data.frame(state=tolower(listeetat[indice]), 
+ party=GV$Party[indice])

On a la visualisation suivante

> library(maps) 
> library(RColorBrewer) 
> couleur=rev(brewer.pal(6, "RdBu")) 
> Z=rep(6,length(basegv$party)) 
> Z[basegv$party=="Democratic"]=1 
> VL0=strsplit(map("state")$names,":")
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,basegv$state)
> niveau=Z[ETAT] 
> map("state", fill = TRUE, col=couleur[niveau])

For those who think more variates should be added to the dataset, some can be found e.g. on http://www.statemaster.com/, like the total executions since 1930, or the date the state joint the U.S.A.

Régression, variables explicatives et géométrie

La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections

Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW15.gif sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW17.gif, deux ensembles (a priori disjoints) de variables explicatives,


Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont


de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par


La première partie correspond à la régression de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW15.gif sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW13.gif, mais il reste un second terme dès lors que https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW17.gif ne sont pas orthogonales. Notons https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW04.gifla matrice de projection (orthogonale) sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif


Alors on peut écrire


En posant

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW07.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW08.gif,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif , i.e. en transformant les variables),


Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents

Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples


Le graphique ci-dessous au cas non orthogonal,


Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,

> chicago=read.table(
+ "http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
> tail(base)
Y    X1 X2     X3
42  4.8 0.152 19 13.323
43 10.4 0.408 25 12.960
44 15.6 0.578 28 11.260
45  7.0 0.114  3 10.080
46  7.1 0.492 23 11.428
47  4.9 0.466 27 13.731
> regression=lm(Y~X1+X2+X3)
> summary(regression)

lm(formula = Y ~ X1 + X2 + X3)

Min     1Q Median     3Q    Max
-9.737 -4.565 -1.479  3.751 16.079

Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X1          -0.62764    5.28130  -0.119 0.905953
X2           0.22378    0.06161   3.632 0.000744 ***
X3          -1.55059    0.38195  -4.060 0.000204 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.527 on 43 degrees of freedom
Multiple R-squared: 0.4417,	Adjusted R-squared: 0.4027
F-statistic: 11.34 on 3 and 43 DF,  p-value: 1.314e-05

> n=length(Y)
> X=matrix(c(rep(1,n),X1,X2,X3),n,4)
> X[1:5,]
[,1]  [,2] [,3]   [,4]
[1,]    1 0.604   29 11.744
[2,]    1 0.765   44  9.323
[3,]    1 0.735   36  9.948
[4,]    1 0.669   37 10.656
[5,]    1 0.814   53  9.730

Ensuite, on va projeter seulement sur les deux premières variables (et la constante)

> FWX1=X[,1:3]
> regression12=lm(Y~X1+X2)
> solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y
[1,]  0.08069764
[2,] 11.56913900
[3,]  0.15108490
> summary(regression12)$coefficients
Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  0.08069764 3.49182154 0.02311047 0.98166664
X1          11.56913900 5.05014782 2.29085156 0.02681947
X2           0.15108490 0.06854408 2.20420043 0.03278652
> FWX2=X[,4]
> H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1)
> M1=diag(rep(1,n))-H1
> FWX2s=M1%*%FWX2
> FWYs =M1%*%Y
> (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs)
[1,] -1.550594
> summary(regression)$coefficients[4]
[1] -1.550594
> (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y-
+        solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2)
[1,] 22.0752495
[2,] -0.6276442
[3,]  0.2237765
> summary(regression)$coefficients[1:3]
[1] 22.0752495 -0.6276442  0.2237765

On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.

Tail index estimation

These data were collected at Copenhagen Reinsurance and comprise 2167 fire losses over the period 1980 to 1990, They have been adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Kron. Note that it is possible to work with the same data as above but the total claim has been divided into a building loss, a loss of contents and a loss of profits.

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> base2=read.table(
+ "http://freakonometrics.free.fr/danish-multivariate.txt",
+ header=TRUE)

Consider here the first dataset (we deal – so far – with univariate extremes),

> X=base1$Loss.in.DKM
> D=as.Date(as.character(base1$Date),"%m/%d/%Y")
> plot(D,X,type="h")

The graph is the following,

A natural idea is then to plot



> Xs=sort(X)
> logXs=rev(log(Xs))
> n=length(X)
> plot(log(Xs),log((n:1)/(n+1)))

Points are on a straight line here. The slope can be obtained using a linear regression,

> B=data.frame(X=log(Xs),Y=log((n:1)/(n+1)))
> reg=lm(Y~X,data=B)
> summary(reg)

lm(formula = Y ~ X, data = B)

Min       1Q   Median       3Q      Max
-0.59999 -0.00777  0.00878  0.02461  0.20309

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.089442   0.001572   56.88   <2e-16 ***
X           -1.382181   0.001477 -935.55   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04928 on 2165 degrees of freedom
Multiple R-squared: 0.9975,	Adjusted R-squared: 0.9975
F-statistic: 8.753e+05 on 1 and 2165 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-500):n,])
> summary(reg)

lm(formula = Y ~ X, data = B[(n - 500):n, ])

Min       1Q   Median       3Q      Max
-0.48502 -0.02148 -0.00900  0.01626  0.35798

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.186188   0.010033   18.56   <2e-16 ***
X           -1.432767   0.005105 -280.68   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07751 on 499 degrees of freedom
Multiple R-squared: 0.9937,	Adjusted R-squared: 0.9937
F-statistic: 7.878e+04 on 1 and 499 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-100):n,])
> summary(reg)

lm(formula = Y ~ X, data = B[(n - 100):n, ])

Min       1Q   Median       3Q      Max
-0.33396 -0.03743  0.02279  0.04754  0.62946

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.67377    0.06777   9.942   <2e-16 ***
X           -1.58536    0.02240 -70.772   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1299 on 99 degrees of freedom
Multiple R-squared: 0.9806,	Adjusted R-squared: 0.9804
F-statistic:  5009 on 1 and 99 DF,  p-value: < 2.2e-16

The slope here is somehow related to the tail index of the distribution. Consider some heavy tailed distribution, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill03.gif, so that https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill27.gif, where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill28.gif is some slowly varying function. Equivalently, the exists a slowly varying function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill29.gif such that https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill30.gif. Then


i.e. since a natural estimator for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill35.gif is the order statistic https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill36.gif, the slope of the straight line is the opposite of tail index https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill98.gif. The estimator of the slope is (considering only the https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill99.gif largest observations)


Hill‘s estimator is based on the assumption that the denominator above is almost 1 (which means that  https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif, as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif), i.e.


Note that, if https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif, but not two fast, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif, then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill12.gif (one can even get https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill11.gif  with stronger convergence assumptions). Further


Based on that (asymptotic) distribution, it is possible to get a (asymptotic) confidence interval for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill98.gif

> xi=1/(1:n)*cumsum(logXs)-logXs
> xise=1.96/sqrt(1:n)*xi
> plot(1:n,xi,type="l",ylim=range(c(xi+xise,xi-xise)),
+ xlab="",ylab="",)
> polygon(c(1:n,n:1),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(1:n,xi+xise,col="red",lwd=1.5)
> lines(1:n,xi-xise,col="red",lwd=1.5)
> lines(1:n,xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to work with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill06.gif, then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill05.gif. And similarly https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill13.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif (and again https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill10.gif with additional assumptions on the rate of convergence), and


(obtained using the delta-method). Again, we can use that result to derive (asymptotic) confidence intervals

> alpha=1/xi
> alphase=1.96/sqrt(1:n)/xi
> YL=c(0,3)
> plot(1:n,alpha,type="l",ylim=YL,xlab="",ylab="",)
> polygon(c(1:n,n:1),c(alpha+alphase,rev(alpha-alphase)),
+ border=NA,col="lightblue")
> lines(1:n,alpha+alphase,col="red",lwd=1.5)
> lines(1:n,alpha-alphase,col="red",lwd=1.5)
> lines(1:n,alpha,lwd=1.5)
> abline(h=0,col="grey")

The Deckers-Einmahl-de Haan estimator is


where for


Then (given again conditions on the speed of convergence i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif, with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif),


Finally, Pickands‘ estimator


it is possible to prove that, as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif,


Here the code is

> Xs=rev(sort(X))
> xi=1/log(2)*log( (Xs[seq(1,length=trunc(n/4),by=1)]-
+ Xs[seq(2,length=trunc(n/4),by=2)])/
+ (Xs[seq(2,length=trunc(n/4),by=2)]-Xs[seq(4,
+ length=trunc(n/4),by=4)]) )
> xise=1.96/sqrt(seq(1,length=trunc(n/4),by=1))*
+sqrt( xi^2*(2^(xi+1)+1)/((2*(2^xi-1)*log(2))^2))
> plot(seq(1,length=trunc(n/4),by=1),xi,type="l",
+ ylim=c(0,3),xlab="",ylab="",)
> polygon(c(seq(1,length=trunc(n/4),by=1),rev(seq(1,
+ length=trunc(n/4),by=1))),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(seq(1,length=trunc(n/4),by=1),
+ xi+xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),
+ xi-xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to use maximum likelihood techniques to fit a GPD distribution over a high threshold.

> library(evd)
> library(evir)
> gpd(X,5)
[1] 2167

[1] 5

[1] 0.8827873

[1] 254

[1] "ml"

xi      beta
0.6320499 3.8074817

xi      beta
0.1117143 0.4637270

[,1]        [,2]
[1,]  0.01248007 -0.03203283
[2,] -0.03203283  0.21504269

[1] "observed"

[1] 0

[1] 754.1115

[1] "gpd"

or equivalently (or almost)

> gpd.fit(X,5)
[1] 5

[1] 254

[1] 0

[1] 754.1115

[1] 3.8078632 0.6315749

[1] 0.1172127

[1] 0.4636270 0.1116136

The interest of the latest function is that it is possible to visualize the profile likelihood of the tail index,

> gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)


> gpd.profxi(gpd.fit(X,20),xlow=0,xup=3)

Hence, it is possible to plot the maximum likelihood estimator of the tail index, as a function of the threshold (including a confidence interval),

> GPDE=Vectorize(function(u){gpd(X,u)$par.ests[1]})
> GPDS=Vectorize(function(u){
+ gpd(X,u)$par.ses[1]})
> u=c(seq(2,10,by=.5),seq(11,25))
> XI=GPDE(u)
> plot(u,XI,ylim=c(0,2))
> segments(u,XI-1.96*XIS,u,XI+
+ 1.96*XIS,lwd=2,col="red")

Finally, it is possible to use block-maxima techniques.

> gev.fit(X)
[1] 0

[1] 3392.418

[1] 1.4833484 0.5930190 0.9168128

[1] 0.01507776 0.01866719 0.03035380

The estimator of the tail index is here the last coefficient, on the right.
Since it is rather difficult to install a package in class rooms, here is the source of rcodes used here (to fit a GPD for exceedances)

> source("http://freakonometrics.blog.free.fr/public/code/gpd.R")

Next time, we will discuss how to use those estimators.