# 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 of an i.i.d. sample  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 =200 observations generated from a Pareto distribution

(here =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 >1.
Since we would like to test  (against the assumption that the expected value is finite, i.e. ), it is natural to consider the likelihood ratio

Under the null hypothesis, the distribution of  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 , and thus, to fit a GPD on the largest 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 (), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under the likelihood ratio statistics has a chi-square distribution). If =2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he largest values (and a confidence interval for the estimator), If =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 gets too close from 1, we should consider thresholds large enough). On the other hand, if =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 =2, we clearly rejct the assumption of infinite mean, and similarly, if =1.5 (finite mean, but infinite variance) Here the test is more robust than the one based on the GPD. And if =1 (i.e. infinite mean), again we accept , Note that if =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, In order to describe how EM algorithm works, assume first that both and are 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 so that we write the log likelihood as which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate (that cannot be observed), taking value when is drawn from , and 0 if it is drawn from . More generally (especially in the case we want to extend our model to 3, 4, … mixtures), and . With that notation, the likelihood becomes and the log likelihood 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 is known, denoted . Then I can predict the value of (i.e. and ) for all observations, So I can inject those values into my log likelihood, i.e. in having maximum (no need to run numerical tools here) that will be denoted . And I can iterate from here. Formally, the first step is where we calculate an expected (E) value, where is the best predictor of given my observations (as well as my belief in ). Then comes a maximization (M) step, where using , I can estimate probability . • A more general framework, all parameters are now unkown So far, it was simple, since we assumed that and were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate . Recall that we had which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from . If , instead of being in the segment was in , then we could have considered mean and standard deviations of observations such that =0, and similarly on the subset of observations such that =1. But we can’t. So what can be done is to consider as the weight we should give to observation i when estimating parameters of , and similarly, 1-would be weights given to observation i when estimating parameters of . So we set, as before and then and for the variance, well, it is a weighted mean again, 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 and then 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 , or . 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 proportional to The true model is , and 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 or , depending whether or 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 ), 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. and 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  sa position à la date , puis poser

autrement dit la probabilité de se retrouver à sa position initiale à la date n (et de retrouver ainsi ses clés), le  pouvait être en n’importe quelle dimension, on posera

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

qui est la probabilité de revenir un jour au point de départ, c’est à dire … 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  issue de  repasse un jour par . S’il va à droite ou à gauche avec probabilités  et  respectivement, notons que , car il faut un nombre de pas pairs pour revenir à l’origine, et

Il est possible de montrer que

et en particulier

On peut relier les  et les  en notant que

Or  constitue une distribution de probabilité, dont on peut définir la fonction génératrice

On notera que que l’on peut écrire

soit, en utilisant l’expression donnée auparavant

Soit finalement

ce qui se traduit par  si la série  converge. Aussi,

donc si , notre alcoolique repassera presque sûrement par le point d’origine, alors que si , il ne repasse pas par le bar où il a oublié ses clés avec probabilité |.

• 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  issue de  repasse un jour par . 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
, car il faut un nombre de pas pairs pour revenir à l’origine, et

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

et que

aussi,

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

La série  est alors divergente, et donc . 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  issue de  repasse un jour par . 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

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

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

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 ?