Detecting distributions with infinite mean

In a post I published a few month ago (in French, here, based on some old paper,there), I mentioned a statistical procedure to test if the underlying distribution https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO02.pngof an i.i.d. sample https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO01.png had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO03.png=200 observations generated from a Pareto distribution

https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO04.png

(here https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, as a start)

> a=2
> X=runif(200)^(-1/a)

Here, we will use the package developped by Mathieu Ribatet,

> library(RFA)
  • Using Generalized Pareto Distribution (and LR test)

A first idea is to fit a GPD distribution on observations that exceed a threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO11.png>1.
Since we would like to test https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png (against the assumption that the expected value is finite, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO07.png), it is natural to consider the likelihood ratio

https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO08.png

Under the null hypothesis, the distribution of https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO09.png should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let  us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill’s estimator) to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO10.png, and thus, to fit a GPD on the https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.pnglargest values. And then to plot all that on a graph (like the Hill plot)

> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }

Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model (https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png the likelihood ratio statistics has a chi-square distribution).
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png largest values (and a confidence interval for the estimator),

If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance), we have

i.e. for those two models, we clearly reject the assumption of infinite mean (even if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png gets too close from 1, we should consider thresholds large enough). On the other hand, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),

  • Using Hill’s estimator

An alternative could be to use Hill’s estimator (with Alexander McNeil’s package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill’s estimator,

> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)

Again, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we clearly rejct the assumption of infinite mean,

and similarly, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance)

Here the test is more robust than the one based on the GPD. And if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean), again we accept https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png,

Note that if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=0.7, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,

(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

Tiens voilà, la pluie. Ah! quel sale temps. Où est-il l’été ? l’été où est-il ?

Dans Les cafards, de Jo Nesbø, Harry Hole (qui vient d’Oslo pour ceux qui n’ont pas dévoré la série) s’entend dire à un moment qu’à Bangkok, les discussions quotidiennes portent rarement sur “la pluie et le beau temps” pour la bonne raison que le climat est assez prévisible à Bangkok (par contre on parle quotidiennement du trafic routier). J’avais évoqué ici l’idée reçue que nous avions, sur le fait que le temps à Rennes n’est pas si changeant que ça. Et je me suis rendu compte depuis que nous sommes à Montréal que le temps pouvait vraiment changer d’un jour sur l’autre (et encore, paraît-il, je n’ai rien vu).
Alors tout d’abord, pour ceux qui en douteraient, la température à Rennes, et à Montréal, ça n’est pas vraiment la même chose,

pour Montréal, alors que pour Rennes on obtient (je précise mais je pense qu’on aurait deviné sans)

Autrement dit, on a une dispersion presque deux fois plus grande à Montréal. Si regarde également les variations (min/max) dans la journée, on obtient, pour Montréal les courbes suivantes,

avec enrouge, la valeur moyenne (sur 15 ans) du maximum observé dans la journée, et enbleu, la valeur moyenne du minimum. A Rennes, on retrouve le fait que la dispersion est relativement faible,

En fait, les choses sont encore pire quand on regarde les quantiles, i.e. le pire minimum observé dans 10% des cas (les plus froids) et le pire maximum observé dans 10% des cas (les plus chauds);

à  Montréal, alors qu’à Rennes, on est beaucoup plus resserrés,

Mais ces pires de cas ne sont pas obtenus dans la même journée, on ne peut pas en dire grand chose. Si on veut aller un peu plus loin, on peut regarder la dispersion dans la journée (i.e. différence entre le minimum et le maximum). On note qu’il est assez stable à Montréal (toujours d’environ 10 degrés)

alors qu’à Rennes, la différence est de 5 degrés l’hiver, contre 10 degré l’été,

Mais surtout, si on regarde la variation d’un jour sur l’autre, danns 90% des cas, on varie de +/- 8 degrés, sur la température moyenne journalière,

avec pas mal de cas un peu extrêmes, avec des variation de 10 ou 15 degré d’un jour sur l’autre sur la température moyenne, alors qu’à Rennes, on a plutôt +/- 5 degrés (avec très peu d’évènements extrêmes, car j’ai tronqué les ordonnées afin de mieux voir au centre)

Bref, je ne suis pas prêt de m’arrêter de parler longtemps de la pluie et du beau temps !

EM and mixture estimation

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM01.png

In order to describe how EM algorithm works, assume first that both https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.pngare perfectly known, and the mixture parameter is the only one we care about.

  • The simple model, with only one parameter that is unknown

Here, the likelihood is

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM04.png

so that we write the log likelihood as

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM05.png

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (that cannot be observed), taking value when https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM07.png is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and 0 if it is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png. More generally (especially in the case we want to extend our model to 3, 4, … mixtures), https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png.
With that notation, the likelihood becomes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM10.png

and the log likelihood

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM11.png

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png is known, denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png. Then I can predict the value of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png) for all observations,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM14.png

So I can inject those values into my log likelihood, i.e. in

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM15.png

having maximum (no need to run numerical tools here)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM16.png

that will be denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM17.png. And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngis the best predictor of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM19.png given my observations (as well as my belief in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png). Then comes a maximization (M) step, where using https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png, I can estimate probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png.

  • A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM30.png. Recall that we had https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png.
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png, instead of being in the segment https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM31.png was in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM32.png, then we could have considered mean and standard deviations of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=0, and similarly on the subset of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=1.
But we can’t. So what can be done is to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png as the weight we should give to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and similarly, 1-https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngwould be weights given to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png.
So we set, as before

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM33.png

and then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM34.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM35.png

and for the variance, well, it is a weighted mean again,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM36.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM37.png

and this is it.

  • Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png and then https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }
Then I get https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png, or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png. So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?… next time, maybe….

Optimization and mixture estimation

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised…

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-01.png proportional to

https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-02.png

The true model is https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-03.png, and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+  sd2 <- param[5]
+  -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 – p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+  }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+   fn = logvrainegX)
Actually, since I might have identification problems, I take either https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-06.png, depending whether https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-07.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-08.png is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: theblack vertical segments. Obviously, when the sample is not enough heterogeneous (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-09.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical  starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again… so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad…. but the probability to be far away from the tru value is not small at all… except when the difference between the two means exceeds 3…
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad… and the two algorithm perform similarly (maybe the first one being a bit closer to the trueparameter).

La tarification avec SAS

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble  que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

DATA contrats;
SET lib.contrats;
lnexpo = log(expo);
run;

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1_Param
           Type3=Genmod1_Var
           Modelfit=Genmod1_InfoModele; 
MODEL nbsin = ageconducteur /
                  dist = poisson   
                  link = log   
                  offset = lnexpo 
                  type3;
RUN; QUIT;

La sortie SAS a alors l’allure suivante

                                  The GENMOD Procedure
                    Critère pour évaluer la qualité de l'ajustement
              Critère                   DF          Valeur       Valeur/DF
              Deviance                63E3      26872.5334          0.4237
              Scaled Deviance         63E3      26872.5334          0.4237
              Pearson Chi-Square      63E3      73275.5362          1.1553
              Scaled Pearson X2       63E3      73275.5362          1.1553
              Log Likelihood                   -18474.2667

       Algorithm converged.
                       Analyse des résultats estimés de paramètres

                              Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

                         Statistiques LR pour Analyse de Type 3
                      Source           DF      Khi 2    Pr > Khi 2
                      ageconducteur     1     148.72        <.0001

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

PROC GAM DATA = base;
MODEL nbsin = spline(ageconducteur) / dist = Poisson;
OUTPUT OUT=gam PREDICTED; 
RUN; QUIT;
PROC SORT DATA = gam NODUPKEY; BY age_cond; RUN; QUIT;

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;
RUN;

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
PROC MEANS DATA = base NOPRINT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo
y_uni_nbsin;
RUN; QUIT;

ensuite, on fait un modèle GLM, et  un modèle GAM, et on récupère les sorties

PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
      gam;
BY age_cond;
RUN;

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1bis_Param
           Type3=Genmod1bis_Var
           Modelfit=Genmod1bis_InfoModele; 
MODEL nbsin = ageconducteur /
                 dist = poisson  
                 link = log      
                 offset = lnexpo 
                 type3           
                 scale = deviance;
RUN; QUIT;

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

                      Analyse des résultats estimés de paramètres

                                   Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

%MACRO CALCUL_AIC_BIC(infomodel=, param=);
    DATA _null_;
    SET &infomodel.;
    IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
    IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
    RUN;
    DATA _null_;
    SET &param.  end=fin;
    RETAIN nb_df 0;
    nb_df = nb_df + df;
    IF fin THEN CALL SYMPUT("k", nb_df);
    RUN;
    DATA genmod_aic_bic;
    SET &param.;
    FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
    Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
    IF Parameter = "Intercept";
    KEEP Loglike K N AIC_CALC BIC_CALC;
    RUN;
    DATA genmod_aic_bic;
    SET genmod_aic_bic;
    Loglike = &loglike.;
    K = &k.;
    N = %eval(&n_etoile. + &k.);
    AIC_CALC = 2 * Loglike + 2 * K;
    BIC_CALC = 2 * Loglike + K * log(N);
    RUN;
    PROC PRINT DATA = genmod_aic_bic;
    RUN; QUIT;
%MEND CALCUL_AIC_BIC;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);

Qui peut m’aider à comprendre les sorties de SAS ?

Je m’étais promis que j’évoquerais une bizarrerie rencontrée avec SAS lors d’une formation…. Écrire ce billet permettra à ceux qui auraient des éléments d’explication de poster un commentaire.
Pour cela, comparons une régression logistique faite avec deux outils différents, sous SAS,

  • avec la procédure logistique

Le code pour faire une régression logistique ressemble à ça

PROC LOGISTIC DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / selection=stepwise;
RUN; QUIT;

ce qui donne la sortie suivante (je passe l’introduction pour insister sur les coefficients)

                                 The LOGISTIC Procedure

                   Analyse des estimations de la vraisemblance maximum
                                                     Erreur         Khi 2
  Paramètre                    DF    Estimation         std       de Wald    Pr > Khi 2

  Intercept                     1        1.7833      0.0676      696.9022        <.0001
  sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
  age_soc      1_AGESOC_-60     1        0.4578      0.0667       47.1020        <.0001
  fract_paiemt Annuel           1        0.6021      0.0997       36.4862        <.0001
  fract_paiemt Mensuel          1       -0.5410      0.0842       41.2342        <.0001
  • avec la procédure genmod (car la régression logistique est un glm)

On peut faire exactement la même chose (théoriquement) en ajustement un modèle GLM,

PROC GENMOD DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / dist = binomial;
RUN;

et la sortie ressemble à ça

                                  The GENMOD Procedure
                      Analyse des résultats estimés de paramètres

                                                  Erreur      Wald 95Limites
Paramètre                     DF   Estimation   standard      de confiance %       Khi 2
Intercept                      1       1.5073     0.1501     1.2131     1.8014    100.85
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000       .
age_soc        1_AGESOC_-60    1       0.9156     0.1334     0.6542     1.1771     47.10
age_soc        Z_AGESOC_+60    0       0.0000     0.0000     0.0000     0.0000       .
fract_paiemt   Annuel          1       0.6634     0.1770     0.3165     1.0104     14.05
fract_paiemt   Mensuel         1      -0.4798     0.1510    -0.7759    -0.1838     10.09
fract_paiemt   Semestriel      0       0.0000     0.0000     0.0000     0.0000       .
Scale                          0       1.0000     0.0000     1.0000     1.0000
  • comparaison des deux sorties

Si on regarde l’impact du sexe par exemple, dans la première sortie on peut lire

sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
alors que dans la seconde sortie, on a
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000

On dira ce qu’on veut, mais moi je trouve cette différence troublante…. Dans la seconde sortie, le coefficient vaut le double de l’autre….
Alors SAS semble s’y retrouver car si on lui demande d’afficher le score prédit pour un individu au hasard (le premier de la base par exemple), les prédictions sont très proches,

                           fract_                                 proba1_       proba1_
  Obs  sexe_soc   age_soc  paiemt      SPOCAM  proba1_logit        
1      Homme          71  Annuel         0      0.10242637    0.10241302

Si quelqu’un sait interpréter ce qui est fait avec cette procédure logistique (car R donne la même chose que la sortie GLM), je suis preneur…..


“Le prix de l’or bat de nouveaux records”… et alors ?

Histoire de changer un peu, un rapide billet d’économie… bien que l’économie monétaire ne soit pas mon domaine de prédilection, loin de là. Mes compétences en la matière (comme sur beaucoup de notions en économie) se limitent à la lecture (mais aussi une intense relecture) d’Obélix et Compagnie

je dois avouer que je suis surpris par les commentaires de certains journalistes économiques, qui nous expliquent que “l’or a atteint des niveaux historiques” (ici, même si ça date un peu) ou que “le prix de l’or bat de nouveaux records” (ce matin cette fois, ). On peut aussi lire ici ou . On peut s’en convaincre en jetant un œil aux graphiques ici, , ou encore. Sans vouloir être méchant, si je transpose ca a une discussion que je pourrais avoir avec mes enfants, ça donnerait: “On a mis 150 euros sur un livret A le jour de votre naissance… c’est incroyable, mais tous les 15 jours, on bat des records: vous êtes de plus en plus riche !”

En fait, effectivement, le prix de l’or atteint des records… mais qu’est-ce que je peux faire de ce “prix historique” quand on sait qu’à coté, tout augmente (comme dirait ma grand mère)… On pourrait normaliser par rapport a autre chose, par exemple un indice boursier (car on parle d’investissement), ou comme le font plutôt les économistes, avec quelque chose en rapport avec la masse monétaire. Non ?
Faisons les calculs: les séries du prix de l’or se trouvent ici par exemple (avec des séries journalières pour les années les plus récentes), et pour les données permettant d’actualiser, on pourra aller ici. Pour faire simple, je me suis contenté des données mensuelles en ligne ici, et pour la masse monétaire, M3 serait pas mal (je renvoie ici pour les détails), mais la Réserve Fédérale américaine ne le publie plus depuis 2006. J’aurais bien aimé aller ici, mais c’est payant. Finalement je suis allé , en extrapolant M3 sur les années récentes (de 2006 à 2010) en utilisant le taux de croissance observé sur M2. Les données ne sont pas corrigées des variations saisonnières.

>  OR = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/prix-or-mensuel.csv",
+ sep=";",header=TRUE)
>  tOR=as.vector(t(OR[,2:13]))
>  VOR=tOR[1:511]

car j’avais mis des 0 pour les mois non observés de 2010 (on les supprime tout simplement). Sinon on retrouve la courbe qui indique qu’effectivement, on a atteint des records…

>    TOR=ts(VOR, frequency = 12, start = c(1968, 1))
>    plot(TOR,col="red")

So what ? Regardons un peu notre masse monétaire qui permettra de ramener le prix de l’or dans le passé au prix réel de l’argent,

>    M = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/M1M2M3.txt")
>    M2=as.numeric(as.character(rev(M$V4)))
>    M3=as.numeric(as.character(rev(M$V5)))
>    M2b=M2[560:length(M2)]/M2[559]
>    M3[560:length(M3)]=M3[559]*M2b
>    TM3=ts(M3,frequency=12,start=c(1959,1))
>    plot(TM3,col="blue")

Bon, j’ai fait une approximation rapide, mais c’est vraiment histoire de voir ce qui se passe si on normalise le prix de l’or. Ceux qui veulent aller plus loin peuvent améliorer le code. Et en corrigeant pour travailler à quantité d’argent équivalent, on obtient des choses un peu différentes (j’aime les euphémismes)

>    TM3b=ts(M3[(12*9+1):length(M3)],
+    frequency = 12, start = c(1968, 1))
>    TORb=TOR/TM3b
>    plot(TORb,col="purple")
>    abline(h=TORb[length(TORb)],lty=2)

On a certes une hausse relativement importante depuis un peu moins de 10 ans du prix de l’or, mais le niveau n’est en rien comparable avec la crise qui s’est observé en 1980 (évoqué ici, , ou encore). Fort heureusement, on ne nous tient pas le même discours sur le PIB, nous expliquant qu’il ne cesse de battre, mois après mois, de nouveaux records ! On notera que ce graphique se rapproche d’autres, en ligne, ici ou
Pour ceux qui ne sont pas convaincu que l’or puisse être relié à une grandeur monétaire, on peut voir l’or comme une alternative à un placement en bourse, et comparer à un indice boursier (par exemple le SP500 pour rester aux États Unis),

>    SP=read.table("http://ichart.finance.yahoo.com/
table.csv?s=%5EGSPC&a=00&b=1&c=1968&d=08
&e=18&f=2010&g=m&ignore=.csv",
+       sep=",",header=TRUE)
>    tSP=rev(SP$Close)
>    TSP=ts(tSP, frequency = 12, start = c(1968, 1))
>    plot(TSP,col="red")

>    TORc=TOR/TSP
>    plot(TORc,col="purple")
>    abline(h=TORc[length(TORc)],lty=2)

Bref, quelle que soit la quantité retenue pour normaliser, on se rend compte que les niveaux atteints ces derniers mois sont très loin de ceux observés en 1980, voire dans les années 80. Pas de quoi en faire la une de la presse….

Marches aléatoires en dimension 1, 2 ou 3, ou histoire éthylique

Revenons à nos histoires d’alcoolique un instant (les précédentes sont en ligne iciet ). L’avant dernière fois (ici), nous avions parlé de la marche aléatoire, et du temps mis avant de trouver la sortie de la rue de la soif.

Et j’avais évoqué la probabilité la probabilité qu’il ne récupère pas les clés (i.e. qu’il ne repasse pas par l’origine) s’il peut sortir un jour de la rue. En probabilité, on appelle cela la probabilité de repasser par l’origine et cela vérifie une propriété étrange dès lors que l’on change de dimension… On va suivre notre alcoolique et noter https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy01.png sa position à la date https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy02.png, puis poser

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy03.png

autrement dit la probabilité de se retrouver à sa position initiale à la date n (et de retrouver ainsi ses clés), le https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png pouvait être en n’importe quelle dimension, on posera

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy05.png

correspondant à la probabilité de premier retour à la position initiale, et enfin

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy06.png

qui est la probabilité de revenir un jour au point de départ, c’est à dire https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png… C’est cette dernière grandeur qui nous intéresse….

  • De la probabilité de récupérer ses clés dans la rue

La rue est ici infinie (contrairement à l’autre fois, où on pouvait en sortir, soit vers la place de Lices, soit vers la place Sainte Anne). Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy08.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. S’il va à droite ou à gauche avec probabilités https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl99.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl98.png respectivement, notons que https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy11.png, car il faut un nombre de pas pairs pour revenir à l’origine, et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy12.png

Il est possible de montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy13.png

et en particulier

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy14.png

On peut relier les https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy15.png et les https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy16.png en notant que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy17.png

Or https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy16.png constitue une distribution de probabilité, dont on peut définir la fonction génératrice

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl18.png

On notera que que l’on peut écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy19.png

soit, en utilisant l’expression donnée auparavant

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl19.png

Soit finalement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl20.png

ce qui se traduit par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy21.png si la série https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy22.png converge. Aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl25.png

donc si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy26.png, notre alcoolique repassera presque sûrement par le point d’origine, alors que si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl97.png, il ne repasse pas par le bar où il a oublié ses clés avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy27.png|.

  • De la probabilité de récupérer ses clés sur une place

La place est ici infinie de part et d’autre. Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy09.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. Il a alors 1 chance sur 4 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite. On peut noter que là encore
https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy11.png, car il faut un nombre de pas pairs pour revenir à l’origine, et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl26.png

Or on peut réécrire cette somme en notant que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy28.png

et que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl29.png

aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl30.png

que l’on peut approcher (via la formule de Stirling) sous le forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy31.png

La série https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy22.png est alors divergente, et donc https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl31.png. Aussi, même sur une place, notre alcoolique reviendra presque sûrement à sa place initiale (et ce une infinité de fois, car une fois qu’il y est revenu et qu’il en repart, il repassera presque sûrement à cette position, etc).

  • De la probabilité de récupérer ses clés dans l’espace

Bon, cette fois notre alcoolique passe en dimension 3…. j’admets, il faut de l’imagination. On s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy10.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. Il a alors 1 chance sur 6 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite, d’aller en haut ou d’aller en bas.

Il y a désormais 6 voisins, et en raisonnant comme avant (le maillage est un peu plus complexe), on obtient que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl32.png

Je laisse les amateurs de combinatoires faire des calculs ou des approximations, mais on peut montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl33.png

La série \sum u_n est alors convergente, et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl34.png

qui vaut un peu plus de 1/3. Autrement dit, on a seulement 1 chance sur 3 de revenir à l’origine et de retrouver ses clés…
Bref, en dimension 1 ou 2, on repasse presque sûrement une infinité de fois à l’origine, mais en dimension 3, on a 2 chances sur 3 de ne jamais y repasser…. J’ai toujours trouvé ça étonnant….

Ramadan, entre le sexe et la faim…

Dans mon précédant billet (ici) je volais un dessin de Martin Vidberg qui parlait de google et de statistiques… Google trend est un outil magique pour noter des tendances et faire des statistiques simples et marrantes. Par exemple, puisque le ramadan touche à sa fin, on peut regarder ce qui motivent les algériens pendant le ramadan. Manifestement, le sexe n’est pas une priorité, que ce soit en anglais,

ou en français,

De manière générale, le mot clé le plus tapé sous google a tendance à se faire tout petit pendant la période du ramadan, alias porno,

En fait, on cherche un peu moins les femmes sur le net,

Par contre, visiblement, on a faim, et on veut bien manger: cuisine explose pendant le ramadan,

On peut aussi regarder les autres pays à forte population musulmane, et observer des tendances similaires, comme le Maroc,

les Émirats Arabes Unis,

ou encore le Pakistan,

(dans ce dernier cas, on observe un fort retour du mot clé sexune fois le ramadan passé… une espèce de manque virtuel). Bref, avec Google, on peut vraiment observer des choses intéressantes… non ?