Tag Archives: simulation

La parité, comme condition d’équilibre stable non dégénéré (simulations)

J’avais mis en ligne hier soit un court billet, expliquant pourquoi la parité ne devait pas uniquement être vu sous un angle éthique, ou moral, mais qu’il convenait de regarder la dynamique d’un système. Et moyennant quelques hypothèses simples (simplistes ?) basées sur des faits stylisées de stigmatisation, on pouvait avoir l’intuition que si le taux de personnes issues de la minorité est trop faible, on converge vers un équilibre dégénéré.

C’est ce qu’on peut visualiser sur la figure suivante

sur une simulation de trajectoire sur 150 pas de temps (150 semaines si on veut, ou 150 jours). On pourrait objecter que c’est juste une simulation, mais en fait, ce comportement quand r_0 est trop faible, on le retrouve assez souvent

avec quelque chose qu’on pourrait interpréter comme un décrochage soudain, où quand r_t passe sous un certain seuil, r_{t+h} plonge vers 0, et (si on reste sur l’interprétation garçon-fille), on converge vers un “Boys’ Club” (ce qu’on appelle en maths, un équilibre dégénéré). Mon modèle est assez simple, et on pourrait le perfectionner (mais l’idée des billets de blog est de faire des choses assez simples). Si on regarde quand on commence proche de 35% pour r_0, on a un ensemble hétérogène de trajectoires, soit vertueuses (qui sont que r_t augmente, pour tendre vers un équilibre paritaire),  soit vicieuses, avec r_t baisse, avec le risque de baisser trop, et de retrouver la spirale infernale précédente.

C’est d’ailleurs très exactement ce qu’on observe si on va plus loin dans les simulations

Une fois ce seuil franchit, r_t plonge vers 0.

Je ne parlerais pas de la partie supérieure, parce que le défaut (ou un des défauts) de mon modèle, c’est que je tiens pas compte du fait que si on passe au dessus de 50%, la “minorité” devient la “majorité”, et donc on pourrait imaginer de manière symétrique le même genre de soucis dans la partie supérieure de la Figure. Au lieu de ces deux tendances, on pourrait imaginer une force de rappel autour de 50% pour maintenir cet état vertueux correspondant à la situation paritaire…

Je peux expliquer un peu comment j’ai obtenu ces figure, en décriant mon modèle dynamique de composition genrée dans un réseau social. On considère un ensemble de n individus connectés au sein d’un graphe aléatoire mixte, évoluant au cours du temps discret t = 0,1,2,\dots. Comme j’avais du mal à gérer les entrées et les sorties, on va constamment suivre n individus, avec à chaque pas de temps, un proportion de femmes est notée r_t \in (0,1) qui évoluera dynamiquement. Le nombre de filles et le nombre de garçons seront ici n_F(t) = \lfloor r_t \, n \rfloor, \quad n_M(t) = n - n_F(t).Les ensembles de nœuds féminins et masculins sont notés F_t et M_t.

Structure du réseau : modèle à blocs stochastiques

Le réseau G_t = (V_t, E_t) est tiré selon un modèle à blocs stochastiques (SBM) à deux blocs (femmes et hommes) avec une matrice de probabilités de connexion P(r_t) qui s’écritP(r) =\begin{pmatrix}p_{FF}(r) & p_{FM}(r) \\p_{FM}(r) & p_{MM}(r) \end{pmatrix}. Je vais définir la fonction logistique (sigmoïde) \sigma(x) = (1 + e^{-x})^{-1},et des fonctions de “rareté” :
\phi_F(r) = \frac{1 - r}{r}, \text{ et } \phi_M(r) = \frac{r}{1 - r}.Les probabilités de lien entre individus dans les groupes s’écrivent alors \begin{aligned}p_{FF}(r) &= \sigma(\theta_0 + \eta + 2\gamma\,\phi_F(r)),\\p_{MM}(r) &= \sigma(\theta_0 + \eta + 2\gamma\,\phi_M(r)),\\p_{FM}(r) &= \sigma(\theta_0 + \delta_{\text{mix}} + \gamma(\phi_F(r) + \phi_M(r))).\end{aligned}Les paramètres \eta et \gamma contrôlent respectivement l’homophilie et la sensibilité à la rareté d’un genre.

Je vais poursuivre la remarque introductive que je faisais auparavant: à chaque étape temporelle, le graphe G_t est resimulé intégralement à partir de la distribution P(r_t). Autrement dit, la structure des liens n’a aucune mémoire, le modèle supposant un mélange social complet à chaque période. Bien sûr, ce choix n’est pas réaliste (dans la réalité, les liens persistent), et donc empêche une vision ABM où on pourrait suivre des trajectoires individuelles… mais ça simplifie grandement la dynamique, en évitant d’avoir à suivre explicitement les entrées, sorties et réajustements de connexions individuels, etc. Mais pour mon billet de blog (qui n’a pas d’autre prétention), ce n’est pas gênant, le modèle représente ainsi une sorte d’état moyen du graph, plutôt qu’un réseau concret avec des individus identifiés.

Exposition locale et propension masculine

Pour chaque nœud masculin i \in M_t, on définit N_t(i) l’ensemble des voisins directs de i, et \deg_t(i)=\#\{N_t(i)\} le degré local (le nombre de voisins directs). On introduit alors s_i(t) = \frac{\#\{|N_t(i) \cap F_t\}}{\deg_t(i)},et z_i(t) = \frac{\#\{N_t(i) \cap M_t\}}{\deg_t(i)} = 1 - s_i(t),
et une mesure de proximité “triangulaire” féminine t_i(t) = \left(\frac{\#\{N_t(i) \cap F_t\}}{\deg_t(i)}\right)^2.Alors, j’ai mis un exposant 2 ici, dans cette formule de “mesure triangulaire féminine” pour avoir une “pondération non linéaire”. Autrement dit, ça permet de renforcer l’effet de la présence multiple de femmes dans le voisinage. L’idée est de représenter l’effet de cohésion féminine locale autour d’un individu masculin. Le but n’est pas de quantifier le fait d’avoir une femme dans le voisinage (déjà capturé par s_i), mais bien d’avoir plusieurs femmes connectées entre elles dans le voisinage, ce va rendre la “présence féminine” visible et normative. Autrement dit : si un homme est connecté à une seule femme, l’effet social reste faible (il peut ignorer la norme) alors que s’il est entouré de plusieurs femmes, qui elles-mêmes sont connectées entre elles, la norme collective est beaucoup plus forte. Il permet de capturer l’idée suivante : plus de femmes proches → plus de triangles F–M–F → plus de pression normative.

Dans ce modèle, la propension individuelle d’un homme à avoir un comportement déplacé (une petite phrase sexiste par exemple) est donnée par \pi_i(t)= \pi_0 \, N(r_t)\, e^{-\alpha\,s_i(t) - \alpha_2\,t_i(t)}\, e^{\zeta\,z_i(t)},N est une fonction normative.

Boucle vertueuse et norme sociale

La norme protectrice (ou “boucle vertueuse”) ne s’active que si la proportion de femmes dépasse un seuil r_{\text{protect}} (il faut qu’il y ait assez de femme). On définit ainsi :
{N}(r) =\begin{cases}\exp\!\big(-\beta_{\text{norm}} \, L(r)\big), & \text{si } r > r_{\text{protect}},\\1, &\text{sinon,}\end{cases}L(r) = \sigma\!\big(\kappa (r - r_\star)\big) est une fonction logistique “raide” (qui tend vers la fonction indicatrice), avec \kappa assez grand et r_\star \approx 0.45 le pivot normatif. Ainsi, lorsque r_t > r_{\text{protect}}, la norme sociale réduit la propension masculine moyenne (à avoir un comportement déplacé).

Incidents et sanctions

Chaque arête (i,j) reliant un homme i et une femme j peut donner lieu à un incident i \to j avec probabilité\mathbb{P}(\text{incident } i \to j) = \lambda \, \pi_i(t),\lambda est un paramètre global de fréquence des comportements déplacés.

On note \mathcal{O}_t l’ensemble des auteurs masculins et\mathcal{V}_t celui des victimes féminines. Quand r_t > r_{\text{protect}}, les auteurs sont sanctionnés avec une probabilité s(r_t) = s_{\max}\, L(r_t), ce qui augmente leur probabilité de sortie.

Sorties et recrutement

Les probabilités de sortie (ou d’abandon) sont données, pour les femmes, parq_j^F(t)= \big(\beta_F+ \tau_{\text{inc}}\,\mathbf{1}\{j\in \mathcal{V}_t\}+ \tau_{\text{iso}}\,\mathbf{1}\{s_j(t) < \theta_{\text{iso}}\}\big)\times B(r_t),et pour les hommes q_i^M(t) = \beta_M + \tau_{\text{sanction}} \, \mathbf{1}\{i \in \mathcal{O}_t \text{ sanctionné}\}.On tire ensuite le nombre de départs féminins X_t^F et masculins X_t^M,X_t^F \sim \mathrm{Binomial}(n_F(t), \bar{q}^F_t),\text{ et }X_t^M \sim \mathrm{Binomial}(n_M(t), \bar{q}^M_t).

Les départs sont compensés par un recrutement de nouveaux membres, avec probabilité p_t^{\text{newF}} d’être une femme :p_t^{\text{newF}}=\begin{cases}r_t, & \text{si } r_t \le r_{\text{protect}},\\[6pt]\mathrm{clamp}\big(r_t + k(0.5 - r_t), 0, 1\big),& \text{si } r_t > r_{\text{protect}}.\end{cases}

Loi d’évolution de la proportion féminine

Après sorties et recrutements, la proportion féminine à l’instantt+1 est r_{t+1} = \frac{n_F(t) - X_t^F + Y_t}{n}= r_t + \frac{Y_t - X_t^F}{n},
Y_t est le nombre de femmes recrutées au temps t. Et en prenant l’espérance conditionnelle, on obtient la dynamique moyenne,\mathbb{E}[r_{t+1}\mid r_t]= r_t + \frac{(X_t^F + X_t^M)\,p_t^{\text{newF}}- \sum_{j\in F_t} q_j^F(t)}{n}.

Interprétation qualitative

Le système présente deux régimes distincts :

  • Sous le seuil r_{\text{protect}} : la boucle vertueuse est inactive, les incidents augmentent, les femmes sortent et on observe l’effondrement dont je parlais en introduction
  • Au-dessus du seuil : la norme protectrice et les sanctions s’activent, les sorties féminines diminuent et on a (en moyenne, et en moyenne seulement comme on le voyait sur les figures introductives) une stabilisation autour de la parité.

Ce modèle reproduit une transition de phase sociale : avec un régime “entre-soi masculin” (boys’ club) et un régime mixte stabilisé par la norme collective. Et dans un monde idéal, on peut imaginer une structure hiérarchique, avec une dizaine de “groupes” (qui représentent une formation par exemple, à l’université), des blocs de blocs, avec la contrainte que si une personne quitte un groupe, forcément elle va dans un autre… et donc un déséquilibre fort dans une formation va avoir des conséquences sur tous les autres (d’autant plus si on met là aussi des matrices de transition entre formations). Donc une formation qui a une parité qui s’effondre, c’est malheureux localement, mais potentiellement dévastateur globalement…

Monte Carlo techniques to create counterfactuals

In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample \{y_1,\cdots,y_n\}, and more generally, in econometrics \{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}. But let’s get back to statistics (without covariates) to illustrate. We assume that observations y_i are realizations of an underlying random variable Y_i. We assume that Y_i are i.id. random variables, with (unkown) distribution F_{\theta}. Consider here some estimator \widehat{\theta} – which is just a function of our sample \widehat{\theta}=h(y_1,\cdots,y_n). So \widehat{\theta} is a real-valued number like . Then, in mathematical statistics, in order to derive properties of the estimator \widehat{\theta}, like a confidence interval, we must define \widehat{\theta}=h(Y_1,\cdots,Y_n), so that now, \widehat{\theta} is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, \widehat{\theta} is .

There are two strategies here. In classical statistics, we use probability theorem, to derive properties of \widehat{\theta} (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, \{y_1,\cdots,y_n\}, and that’s a pity. But maybe we can create another one \{y_1^{(1)},\cdots,y_n^{(1)}\}, as realizations of F_{\theta}, and another one \{y_1^{(2)},\cdots,y_n^{(2)}\}, anoter one \{y_1^{(3)},\cdots,y_n^{(3)}\}, etc. From those counterfactuals, we can now get a collection of estimators, \widehat{\theta}^{(1)},\widehat{\theta}^{(2)}, \widehat{\theta}^{(3)}, etc. Instead of using mathematical tricks to calculate \mathbb{E}(\widehat{\theta}), compute \frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}That’s what we’ve seen last friday.

I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students,

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
X=Davis$height[Davis$sex=="M"]

We can visualize its distribution (density and cumulative distribution)

u=seq(155,205,by=.5)
par(mfrow=c(1,2))
hist(X,col=rgb(0,0,1,.3))
lines(density(X),col="blue",lwd=2)
lines(u,dnorm(u,178,6.5),col="black")
Xs=sort(X)
n=length(X)
p=(1:n)/(n+1)
plot(Xs,p,type="s",col="blue")
lines(u,pnorm(u,178,6.5),col="black")

Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a \mathcal{N}(178,6.5^2). In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals

hist(X,col=rgb(0,0,1,.3))
lines(density(X),col="blue",lwd=2)
  Y=rnorm(n,178,6.5)
  hist(Y,col=rgb(1,0,0,.3))
  lines(density(Y),col="red",lwd=2)
Ys=sort(Y)
plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205))
polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA)
lines(Xs,p,type="s",col="blue",lwd=2)
lines(Ys,p,type="s",col="red",lwd=2)

We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution \widehat{F} obtained from \{y_1,\cdots,y_n\} (the blue curve), and some conterfactual, \widehat{F}^{(s)} obtained from \{y_1^{(s)},\cdots,y_n^{(s)}\} generated from F_{\theta_0} – where \theta_0 is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance.

d=rep(NA,1e5)
for(s in 1:1e5){
d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic
}
ds=density(d)
plot(ds,xlab="",ylab="")
dks=ks.test(X,"pnorm",178,6.5)$statistic
id=which(ds$x&gt;dks)
polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA)
abline(v=dks,col="red")

If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic \widehat{d}^{(1)}, \widehat{d}^{(2)}, etc, and compare it with the one observe on our sample \widehat{d}. The proportion of samples where the test-statistics exceeded the one observed

mean(d&gt;dks)
[1] 0.78248

is the computational version of the p-value

ks.test(X,"pnorm",178,6.5)
 
	One-sample Kolmogorov-Smirnov test
 
data:  X
D = 0.068182, p-value = 0.8079
alternative hypothesis: two-sided

I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations).

Graduate Course on Advanced Tools for Econometrics (1)

This Monday, I will be giving the first part of the (crash) graduate course on advanced tools for econometrics. It will take place in Rennes, IMAPP room, and I have been told that there will be a visio with Nantes and Angers. Slides for the morning are online, as well as slides for the afternoon.

In the morning, we will talk about smoothing techniques, and in the afternoon, it will be on simulations and bootstrap techniques.

Tests non-paramétriques et simulations

Lors du dernier cours de statistique, nous avons présenter les tests d’ajustment de lois. Nous avions illustré ces tests à partir de la taille d’individus (déjà utilisé pour présenter l’ajustement de lois, et l’estimation de densité), correspondant à https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_n\}.

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)]
> X=Davis$height
> n=length(X)

On notera https://latex.codecogs.com/gif.latex?%20(x_{i:n}) la statistique d’ordre, au sens où

https://latex.codecogs.com/gif.latex?%20x_{1:n}\leq%20x_{2:n}\leq\cdots\leq%20x_{n-1:n}\leq%20x_{n:n}

Parmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant

> PP=function(Y,F=pnorm){
+   n=length(Y)
+   x=F(sort(Y))
+   y=seq(1/n/2,1-1/n/2,by=1/n)
+   return(data.frame(x=x,y=y))
+ }

qui représente (à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?\left\{F_0(x_{i:n}),\frac{i}{n}\right\}

et celui pour le QQ plot serait

> QQ=function(Y,Q=qnorm){
+   n=length(Y)
+   x=Q(seq(1/n/2,1-1/n/2,by=1/n))
+   y=sort(Y)
+   return(data.frame(x=x,y=y))
+ }

qui représente (toujours à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?%20\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}

où https://latex.codecogs.com/gif.latex?%20F_0 est la loi que l’on cherche à tester, au sens où on a https://latex.codecogs.com/gif.latex?%20H_0:F=F_0 avec comme hypothèse alternative https://latex.codecogs.com/gif.latex?%20H_1:F\neq%20F_0.

Continue reading Tests non-paramétriques et simulations

Statistical Tests: Asymptotic, Exact, ou based on Simulations?

This morning, in our mathematical statistics course, we’ve been discussing the ‘proportion test‘, i.e. given a sample of Bernoulli trials, with , we want to test

against 

A natural test (which can be related to the maximum likelihood ratio test) is  based on the statistic

The test function is here

To get the bounds of the acceptance region, we need the distribution of , under . Consider here a numerical application

n=20
p=.5
set.seed(1)
echantillon=sample(0:1,size=n,
            prob=c(1-p,p),
            replace=TRUE)
  • the asymptotic distribution

The first (and standard idea) is to use the central limit theorem, since

So, under ,

Then  while . The acceptance region is then between the two red lines, below,

T=sqrt(n)*(mean(echantillon)-.5)/
  sqrt(mean(echantillon)*
  (1-mean(echantillon)))
u=seq(-3,3,by=.01)
v=dnorm(u)
plot(u,v,type="l",lwd=2)
abline(v=qnorm(.025),col="red")
abline(v=qnorm(.975),col="red")
abline(v=T,col="blue")

  • the exact distribution

Here we use the fact that

Using transformation of the ‘density’, we can (at least numerically) compute the (exact) distribution of

 

u=seq(-3,3,by=.01)
v=sqrt(.5*(1-.5))*n*dbinom(round(
  (sqrt(.5*(1-.5))*u/sqrt(n)+.5)*n),
  size=n,prob=.5)/sqrt(n)

Here I used a round value, it guess it would be better with a floor function, but here the graph looks symmetric (which is something I like)

abline(v=sqrt(n)*(qbinom(.025,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
abline(v=sqrt(n)*(qbinom(.975,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
lines(u,v,type="s")

  • distribution based on Monte Carlo simulations

Probably more interesting, here we do not use the fact that we might know the distribution of the mean. We just generate random samples, under , and then compute ,

T=rep(NA,1000)
for(i in 1:1000){
x=sample(0:1,size=n,
         prob=c(1-.5,.5),
         replace=TRUE)
m=mean(x)
T[i]=(m-.5)/sqrt(m*(1-m))*sqrt(n)}
lines(density(T),lwd=2)
abline(v=quantile(T,.025),col="red")
abline(v=quantile(T,.975),col="red")

L’espérance est un opérateur linéaire, so what ?

Oui, “l’espérance est un opérateur linéaire”. On n’arrête pas d’insister sur cette propriété dans la plupart des cours de probabilité. En fait, je pense que cette propriété a deux implications importantes. La première est que l’on interprète souvent cette phrase en disant que si http://freakonometrics.blog.free.fr/public/perso5/espl01.gif intégrable, alors pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gifhttp://freakonometrics.blog.free.fr/public/perso5/espl7.gif. La conséquence de cette propriété (qui est juste) est que http://freakonometrics.blog.free.fr/public/perso5/espl09.gif dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).

Mais cette linéarité dit davantage. Elle dit que pour toutes variables aléatoires http://freakonometrics.blog.free.fr/public/perso5/espl01.gifet http://freakonometrics.blog.free.fr/public/perso5/espl2.gif (définies sur le même espace), http://freakonometrics.blog.free.fr/public/perso5/espl3.gif, pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gif. Et ça c’est fort. En particulier, on ne dit pas qu’il faut que http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif.

Pourtant @Arnaud (qui passe souvent à l’occasion sur le blog) m’a laissé un couple de questions par courriel: “la dépendance modifie-t-elle la moyenne ?[…] sur des exemples simples en R, la moyenne est toujours perturbée”. L’idée était que les espérances étaient approchées numériquement par simulation. Or les méthodes de Monte Carlo reposent sur la loi des grands nombres et sur le théorème central limite. En particulier, si http://freakonometrics.blog.free.fr/public/perso5/espl34.gif, et que l’on simule les couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif, alors

http://freakonometrics.blog.free.fr/public/perso5/espl31.gif

autrement dit

http://freakonometrics.blog.free.fr/public/perso5/espl32.gif

i.e. la vitesse de convergence est affectée par la dépendance, pas la valeur vers laquelle on va tendre (car la moyenne tend vers l’espérance, qui est linéaire). En particulier, la convergence sera d’autant plus lente que http://freakonometrics.blog.free.fr/public/perso5/espl33.gif sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.

Par exemple, si on s’amuse à simuler 100 couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif Gaussiens, avec

http://freakonometrics.blog.free.fr/public/perso5/espl11.gif

on a les valeurs suivantes pour la moyenne de ce couple (obtenu en générant 10,000 échantillons de 100 paires), avec en abscisse la corrélation entre les deuxvariables

le trait rouge est la moyenne empirique des moyennes obtenues, et les régions sont délimitées par les quantiles à 5% (et 95%), 10% (et 90%) et 25% (et 75%). En moyenne, les simulations donnent la même chose, peu importe la dépendance. Mais plus la dépendance est forte, plus la variance sera grande, et plus la convergence sera lente (et les simulations auront d’autant plus tendance à s’éloigner – ou se disperser autour – de la vraie valeur). Mais cette réponse est partielle, car on a supposé que la variance était finie. Or cette hypothèse n’est pas nécessaire pour assurer la convergence de la moyenne vers l’espérance (et donc pour utiliser des méthodes de simulations). Que se passerait-il avec des lois de variance infinie ?

Pour autre chose que les vecteurs Gaussiens (par exemple des variables de variance infinie, comme des lois de Pareto) on peut utiliser le code suivant (on continue pour l’instant à utiliser une copule Gaussienne),

library(mnormt)
library(copula)
library(evir)
set.seed(1)
ns=100
R=seq(-.95,.95,by=.05)
M=BINF1=BSUP1=rep(NA,length(R))
for(i in 1:length(R)){
r=R[i]
norm.cop = normalCopula(r, dim = 2)
S=rep(NA,20000)
for(j in 1:20000){
U=rcopula(norm.cop,ns)
X=cbind(qgpd(U[,1], .7 ,1),qgpd(U[,2], .5 ,1))
S[j]=mean(X[,1]+X[,2])}
BINF1[i]=quantile(S,.05)
BSUP1[i]=quantile(S,.95)
M[i]=mean(S)
}
plot(R,M,col="white",ylim=range(c(BINF1,BSUP1)))
polygon(c(R,rev(R)),c(BINF1,rev(BSUP1)),
col="light yellow",border=NA)
lines(R,BINF1,lty=2)
lines(R,BSUP1,lty=2)
lines(R,M,lwd=2,col="red")
esperance=1+1/(1-.7)+1+1/(1-.5)
segments(min(R),esperance,max(R),esperance,col="blue")

Sur cette somme de deux lois de Pareto (comme dans le code ci-dessus), on obtient numériquement

Autrement dit, ici, l’effet dépendance semble avoir un impact plus faible sur l’approximation de la moyenne. On peut aussi essayer de changer la copule, par exemple mettre une copule de type Clayton (avec en abscisse cette fois le tau de Kendall),

(qui présente de la dépendance dans la queue inférieure) ou la copule duale (qui présente de la dépendance dans la queue supérieure)

Bref, la structure de la dépendance ne semble pas trop changer l’estimation de la moyenne (ou plutôt la vitesse de la convergence). On peut aussi se demander s’il y a un effet dimension. Par exemple, si au lieu de simuler des couples, on simule des vecteurs de plus grande taille. Par exemple, un vecteur échangeable en dimension 10, avec des marges normales centrées réduites,

norm.cop =normalCopula(r, dim = 10,dispstr="ex")

Moralité ? effectivement, la dépendance a un impact sur la vitesse de convergence quand on travaille sur des vecteurs aléatoires. Mais de manière assez surprenante, c’est surtout le cas pour la loi normale. Et étrangement, cet impact ne semble pas lié à de la dépendance extrême, ou à la dimension…

Quels résidus et quelle loi simuler ?

Suite à une question sur les résidus, je vais essayer de prendre deux minutes pour reprendre un point que je n’avais pas abordé en cours, faute de temps. Les résidus les plus intéressants en provisionnement (et dans les GLM) sont les résidus de Pearson,

http://freakonometrics.blog.free.fr/public/perso4/res-boot-01.gif

Mais il existe également des résidus dits ajustés, permettant de prendre en compte le fait que le nombre de paramètres est ici relativement grand, ramené au nombre d’observations. On pose alors

http://freakonometrics.blog.free.fr/public/perso4/res-boot-2.gif

A première vue, utiliser l’un ou l’autre devrait donner la même chose si l’on génère des pseudo-triangles, puisque dans un cas, on utiliserait

http://freakonometrics.blog.free.fr/public/perso4/res-boot.gif

et dans l’autre

http://freakonometrics.blog.free.fr/public/perso4/res-boot-3.gif

Bref, multiplier par une constante pour ensuite diviser par la même constante, ça revient au même. Sauf qu’il peut sembler légitime de poser

http://freakonometrics.blog.free.fr/public/perso4/res-boot-05.gif

(on reviendra là dessus tout à l’heure). En attendant, voilà nos deux résidus

> (E=residuals(regp,"pearson"))
1             2             3             4
9.488238e-01  2.404895e-02  1.168421e-01 -1.082940e+00
5             6             7             8
1.302749e-01 -1.007348e-13 -1.128013e+00  2.773332e-01
9            10            11            13
5.669707e-02  8.919633e-01 -2.110748e-01 -1.533031e+00
14            15            16            19
-2.213449e+00 -1.024162e+00  4.237393e+00 -4.899687e-01
20            21            25            26
7.929194e-01 -2.972380e-01 -4.275912e-01  4.140426e-01
31
-6.202125e-15
> n=sum(is.na(Y)==FALSE)
> k=ncol(PAID)+nrow(PAID)-1
> (R=residuals(regp,"pearson")*sqrt(n/(n-k)))
1             2             3             4
1.374976e+00  3.485024e-02  1.693203e-01 -1.569329e+00
5             6             7             8
1.887862e-01 -1.459787e-13 -1.634646e+00  4.018940e-01
9            10            11            13
8.216186e-02  1.292578e+00 -3.058764e-01 -2.221573e+00
14            15            16            19
-3.207593e+00 -1.484151e+00  6.140566e+00 -7.100321e-01
20            21            25            26
1.149049e+00 -4.307387e-01 -6.196386e-01  6.000048e-01
31
-8.987734e-15

L’autre point sur lequel je voulais revenir sur les simulations, est que la loi de Poisson n’est peut-être pas adaptée pour simuler des scénarios de paiements futurs. En effet, si on fait une régression quasipoisson, on voit que le paramètre de surdispersion n’est pas négligeable,

> regqp=glm(Y~as.factor(D)+as.factor(A),
+     data=base,family=quasipoisson(link="log"))
> summary(regqp)

Call:
glm(formula = Y ~ as.factor(D) + as.factor(A), 
family = quasipoisson(link = "log"),
data = base)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3426  -0.4996   0.0000   0.2770   3.9355

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)       8.05697    0.02769 290.995  < 2e-16 ***
as.factor(D)2    -0.96513    0.02427 -39.772 2.41e-12 ***
as.factor(D)3    -4.14853    0.11805 -35.142 8.26e-12 ***
as.factor(D)4    -5.10499    0.22548 -22.641 6.36e-10 ***
as.factor(D)5    -5.94962    0.43338 -13.728 8.17e-08 ***
as.factor(D)6    -5.01244    0.39050 -12.836 1.55e-07 ***
as.factor(A)2002  0.06440    0.03731   1.726 0.115054
as.factor(A)2003  0.20242    0.03615   5.599 0.000228 ***
as.factor(A)2004  0.31175    0.03535   8.820 4.96e-06 ***
as.factor(A)2005  0.44407    0.03451  12.869 1.51e-07 ***
as.factor(A)2006  0.50271    0.03711  13.546 9.28e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Dispersion parameter for quasipoisson family taken to be 3.18623

Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

On peut alors envisager de simuler non pas des lois de Poisson, mais des lois quasipoissons, en reprenant un vieux code,

> rqpois = function(n, lambda, phi, roundvalue = TRUE) {
+ b = phi
+  a = lambda/phi
+ r = rgamma(n, shape = a, scale = b)
+  if(roundvalue){r=round(r)}
+ return(r)
+ }

A partir de ces deux remarques, on peut reprendre le code qui permettait de générer des scénarios de paiements pour les années futures,

> Yp=predict(regp,type="response",newdata=base)
> Rs=rep(NA,20000)
> for(s in 1:20000){
+ serreur=sample(erreur,
+ size=36,replace=TRUE)
+ E=matrix(serreur,6,6)
+ sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6))
+ sbase=data.frame(sY=as.vector(sY),D,A)
+ sbase$sY[is.na(Y)==TRUE]=NA
+ sreg=glm(sY~as.factor(D)+as.factor(A),
+ data=sbase,family=poisson(link="log"))
+ sYp=predict(sreg,type="response",
+ newdata=sbase)
+ sYpscenario=rqpois(36,sYp,phi=3.18623)
+ Rs[s]=sum(sYpscenario[is.na(Y)==TRUE])
+ }

si on regarde la densité de la distribution des paiements futurs, on obtient,

plot(density(Rs))

avec en bleu la distribution obtenue sur les résidus de Pearson brut, et en simulant une loi de Poisson, et en rouge la distribution des provisions avec les deux modifications. On observe clairement que les quantiles (par exemple) ont fortement changé. On peut le voir encore plus précisément sur un box-plot

Si on revient cinq minutes sur ce qu’on vient de faire. Dans le premier rappelons que l’on a besoin d’utiliser, pour prédire les paiements futurs, que

http://freakonometrics.blog.free.fr/public/perso4/siiiim-01.gif

(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.

http://freakonometrics.blog.free.fr/public/perso4/siiiiim-03.gif

Comme on utilise un modèle de Poisson, on suppose aussi que

http://freakonometrics.blog.free.fr/public/perso4/siiiim-5.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso4/siiiiim-08.gif

c’est à dire que

http://freakonometrics.blog.free.fr/public/perso4/siiiim-04.gif

Il faut avoir des résidus centrés et de variance unitaire. A première vue, c’est le cas pour nos résidus de Pearson (à peu près),

> mean(residuals(regp,"pearson"))
[1] -0.02462518
> sd(residuals(regp,"pearson"))
[1] 1.261934

sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur http://freakonometrics.blog.free.fr/public/perso4/siiiim-11.gif (classique en statistique quand on possède http://freakonometrics.blog.free.fr/public/perso4/siiim-10.gif observations)

> sqrt(sum((residuals(regp,"pearson")-
+ mean(residuals(regp,"pearson")))^2)/
+ (length(residuals(regp,"pearson"))-1))
[1] 1.261934

Sauf qu’ici on ne perd pas un degré de liberté (remplacer l’espérance des observations par leur moyenne empirique): en régression, il faut corriger par le nombre de variables explicatives. Et donc, pour avoir un estimateur sans biais de la variance de nos résidus, il faut corriger par un facteur http://freakonometrics.blog.free.fr/public/perso4/siiiim-12.gif. D’où l’utilisation des résidus dits ajustés afin d’avoir des résidus de variance (vraiment) unitaire.
Pour le second point, on utilise le fait qu’en pratique, la dispersion des paiements est plus grande que ce que donnerait un modèle de Poisson. Et il convient d’en tenir compte dans nos simulations. En l’occurrence, on veut que  https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp001.png, mais aussi quehttps://perso.univ-rennes1.fr/arthur.charpentier/latex/qp002.png. On a vu que le paramètre de surdispersion ne pouvait être supposé unitaire. On a alors le choix, entre simuler une loi Gamma de paramètres

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp003.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp004.png

ou bien simuler une loi binomiale négative, de moyenne https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp008.png, et de paramètre de surdispersion https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp009.png telle que la variance s’écrive

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

Aussi, on prend https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp011.png alors que

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

Avec ces méthodes, on obtient quelque chose… qui est programmé sous R,

> BootChainLadder(PAID,10000,"od.pois")
BootChainLadder(Triangle = PAID, R = 10000, 
process.distr = "od.pois")

Latest Mean Ultimate Mean IBNR SD IBNR IBNR 75% IBNR 95%
1  4,456         4,456       0.0     0.0        0        0
2  4,730         4,752      22.2    12.2       28       44
3  5,420         5,456      35.5    15.2       43       64
4  6,020         6,086      65.7    19.8       78      102
5  6,794         6,946     151.9    28.7      170      202
6  5,217         7,364   2,147.4   111.0    2,218    2,339

Totals
Latest:         32,637
Mean Ultimate:  35,060
Mean IBNR:       2,423
SD IBNR:           132
Total IBNR 75%:  2,506
Total IBNR 95%:  2,653

> quantile(Rs,c(.75,.95))
75%  95%
2509 2653

et qui donne les mêmes résultats que ce qu’on vient de reprogrammer…