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

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 à $\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 $(x_{i:n})$ la statistique d’ordre, au sens où $x_{1:n}\leq x_{2:n}\leq\cdots\leq x_{n-1:n}\leq x_{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 $\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 $\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}$ où $F_0$ est la loi que l’on cherche à tester, au sens où on a $H_0:F=F_0$ avec comme hypothèse alternative $H_1:F\neq F_0$. # 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$\boldsymbol{x}=\{x_1,\cdots,x_n\}$, with $X_i\sim\mathcal{B}(p)$, we want to test $H_0:p=p_0$against $H_1:p\neq p_0$ A natural test (which can be related to the maximum likelihood ratio test) is based on the statistic $T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$ The test function is here $\psi(\boldsymbol{x})=\boldsymbol{1}(T(\boldsymbol{x})\notin [c_{1,\alpha},c_{2,\alpha}])$To get the bounds of the acceptance region, we need the distribution of $T(\boldsymbol{X})$, under $H_0$. 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 $\sqrt n\frac{\hat{p} - p}{\sqrt{p (1-p)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$ So, under $H_0$, $\sqrt n\frac{\hat{p} - p_0}{\sqrt{p_0 (1-p_0)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$ Then $c_{1,\alpha}=\Phi^{-1}(\alpha/2)$ while $c_{2,\alpha}=\Phi^{-1}(1-\alpha/2)$. 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 $\sum_{i=1}^n X_i \sim \mathcal{B}(n,p)$ Using transformation of the ‘density’, we can (at least numerically) compute the (exact) distribution of $T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$ 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 $H_0$, and then compute $T(\boldsymbol{X})$, 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 intégrable, alors pour tout et . La conséquence de cette propriété (qui est juste) est que 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 et (définies sur le même espace), , pour tout et . Et ça c’est fort. En particulier, on ne dit pas qu’il faut que et sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre et . 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 , et que l’on simule les couples , alors autrement dit 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 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 Gaussiens, avec 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, 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 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 et dans l’autre 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 (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),
+ 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

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

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

i.e.

c’est à dire que

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  (classique en statistique quand on possède  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 . 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  , mais aussi que. 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

et

ou bien simuler une loi binomiale négative, de moyenne , et de paramètre de surdispersion  telle que la variance s’écrive

Aussi, on prend  alors que

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…