Category Archives: Research

PhD defense on copulas

This Wednesday I will be at Université Paris 1 Sorbonne as a member of the jury of the PhD thesis of Pierre-André Maugis, on conditional correlation and vine copula.

Vine copulas were born in 2002 with thepaper of Tim Bedford and Roger M. CookeVines–a new graphical model for dependent random variables. The idea is to use the following decomposition for a multivariate density

(from Bayes formula, with synthetic notations). Then using the relationship between a bivariate density and its copula (density)

thus

Using again Bayes formula,

and we can write

Since  and , the previous expression becomes

or to stress on the most important part (as I see it)

It is common then to assume that this conditional copula does not depend on the conditioning parameter. The more detailed expression of that joint trivariate density is

The (parametric) inference algorithm is defined in Cooke, Joe and Aas (2010) as follows

The important assumption in vine copula models is that conditional copulas are constant. And this assumption might be relevant in some cases. For instance, in the Gaussian case (the observations have a Gaussian joint distribution – or at least copula – and we fit a vine model with Gaussian bivariate copulas).

The code to fit a vine copula is the following,

> library(CDVine)
> library(mnormt)
> SIGMA=matrix(c(1,.6,.7,.6,1,.8,.7,.8,1),3,3)
> X=rmnorm(n=100000,varcov=SIGMA)
> CDVineSeqEst(dat=X, family = c(1,1,1),
+ type = 1, method = "mle")
$par
[1] 0.6001505 0.7023699 0.6698215
 
$par2
[1] 0 0 0

Note that it is consistent with the following algorithm where conditional copulas are fitted. In the following, for all values of the given component, we wit a Gaussian copula for the conditional remaining pair,

> U=pnorm(X)
> U1U2=U[,1:2]
> U1U3=U[,c(1,3)]
> GaussCop = normalCopula(param=.5, dim = 2)
> U1U2=U[,1:2]
> U1U3=U[,c(1,3)]
> fit12.mpl = fitCopula(GaussCop, U1U2, method="mpl")@estimate
> fit13.mpl = fitCopula(GaussCop, U1U3, method="mpl")@estimate
> fit12.mpl
[1] 0.5984932
> fit13.mpl
[1] 0.7005185
> fit23a=fit23b=rep(NA,99)
> for(i in 4:96){
+ x=i/100
+ C12=pcopula(normalCopula(param=fit12.mpl, dim = 2),U1U2)
+ C13=pcopula(normalCopula(param=fit13.mpl, dim = 2),U1U3)
+ U12=rank(C12)/(nrow(U)+1)
+ U13=rank(C13)/(nrow(U)+1)
+ U23=cbind(U12[abs(U[,1]-x)<.02],U13[abs(U[,1]-x)<.02])
+ V23=cbind(rank(U23[,1])/(nrow(U23)+1),
+ rank(U23[,2])/(nrow(U23)+1))
+ fit23.mpl = fitCopula(GaussCop, V23, method="mpl")@estimate
+ fit23a[i]=fit23.mpl
+ }
> plot(X,fit23a,col="red")

It looks like assuming the conditional copula as constant was a valid assumption here

But note that if the true distribution is not Gaussian, then assuming the conditional copula as constant is not valid anymore (here a trivariate Clayton copula was generated)

Climate change and insurance

I will be in Lyon next Monday to give a talk on “Modeling heat-waves: return period for non-stationary extremes” in a workshop entitled “Changement climatique et gestion des risques“. An interesting reference might be some pages from Le Monde (2010). The talk will be more a discussion about modeling series of temperatures (daily temperatures). A starting point might be the IPCC Third Assessment graph which illustrates the effect on extreme temperatures when (a) the mean temperature increases, (b) the variance increases, and (c) when both the mean and variance increase for a normal distribution of temperature.

I will add here some code used to generate some graphs I will comment. The graph below it the daily minimum temperature,

TEMP=read.table("http://freakonometrics.blog.free.fr/
public/data/TN_STAID000038.txt",header=TRUE,sep=",")
D=as.Date(as.character(TEMP$DATE),"%Y%m%d")
T=TEMP$TN/10
day=as.POSIXlt(D)$yday+1
an=trunc(TEMP$DATE/10000)
plot(D,T,col="light blue",xlab="Minimum
daily temperature in Paris",ylab="",cex=.5)
abline(R,lwd=2,col="red")

We can clearly see an increasing linear trend. But we do not care (too much) here about the increase of the average temperature, but more dispersion, and tails. Here are decenal box-plots

or quantile-regressions

library(quantreg)
PENTESTD=PENTE=rep(NA,99)
for(i in 1:99){
R=rq(T~D,tau=i/100)
PENTE[i]=R$coefficients[2]
PENTESTD[i]=summary(R)$coefficients[2,2]
}
m=lm(T~D)$coefficients[2]
plot((1:99)/100,(PENTE/m-1)*100,type="b")
segments((1:99)/100,((PENTE-2*PENTESTD)/m-1)*100,
(1:99)/100,((PENTE+2*PENTESTD)/m-1)*100,
col="light blue",lwd=3)
points((1:99)/100,(PENTE/m-1)*100,type="b")
abline(h=0,lty=2,col="red")

In order to get a better understanding of the graph above, here are slopes of quantile regressions associated to different probabilities,

The annualized maxima (of minimum temperature, i.e. warmest night of the year)

i.e. the regression of yearly maximas.

tail index of a Generalized Pareto distribution

Instead of looking at observation over a century (the trend is obviously linear), we can focus on seaonal behavior,

B=data.frame(Y=rep(T,3),X=c(day,day-365,day+365),
A=rep(an,3))
library(quantreg)
library(splines)
Q50=rq(Y~bs(X,10),data=B,tau=.5)
Q95=rq(Y~bs(X,10),data=B,tau=.95)
Q05=rq(Y~bs(X,10),data=B,tau=.05)
YP95=predict(Q95,newdata=data.frame(X=1:366))
YP05=predict(Q05,newdata=data.frame(X=1:366))
I=(T>predict(Q95))|(T<predict(Q05))
YP50=predict(Q50,newdata=data.frame(X=1:366))
plot(day[I],T[I],col="light blue",cex=.5)
lines(1:365,YP95[1:365],col="blue")
lines(1:365,YP05[1:365],col="blue")
lines(1:365,YP50[1:365],col="blue",lwd=3)

with on red series from 1900 till 1920, and on purple from 1990 till 2010. If we remove the linear trend, and the seasonal cycle, here are the residuals, assume to be stationary,

on during the year

Obviously, something has been missed,

The graph below is the volatility of the residual series, within the year,

Instead of looking at volatility, we can focus on tails, with tail index per month,

mois=as.POSIXlt(D)$mon+1
Pmax=Dmax=matrix(NA,12,2)
for(s in 1:12){
X=T3[mois==s]
FIT=gpd(X,5)
Pmax[s,1:2]=FIT$par.ests
Dmax[s,1:2]=FIT$par.ses
}
plot(1:12,Pmax[,1],type="b",col="blue",
ylim=c(-.6,0))
segments(1:12,Pmax[,1]+2*Dmax[,1],1:12,Pmax[,1]-
2*Dmax[,1],col="light blue",lwd=2)
points(1:12,Pmax[,1],col="blue")
text(1:12,rep(-.5,12),c("JAN","FEV","MARS",
"AVR","MAI","JUIN","JUIL","AOUT","SEPT",
"OCT","NOV","DEV"),cex=.7)

At the end of the talk, I will also mention multiple city models, e.g. Paris and Marseille,

If we look at residuals (once we have removed the linear trend and the seasonal cycle) we observe some positive dependence

In order to study (strong) tail dependence, define

http://freakonometrics.hypotheses.org/files/2017/07/Llatex2png.2.php_.png

for lower left tail and

http://freakonometrics.hypotheses.org/files/2017/07/Clatex2png.2.php_.png

for upper right tail, where http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-12.2.php_.png is the survival copula associated to http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-13.2.php_.png, i.e.
http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-14.2.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-15.2.php_.png

It looks like there is no tail dependence (in the uper tail). But it is also possible to study weaker tail dependence, through

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.3.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.4.php_.png


Slides can be visualized below, I will upload them soon,

Les doctorants, le sexe et le directeur de thèse

Pourquoi certains chercheurs publient beaucoup, et d’autres moins ? ou de l’hétérogénéité dans une loi de Poisson (oui, le but principal du billet est de parler de la régression de Poisson).

Il y a une vingtaine d’années, Scott Long publiait un article intitulé ‘The Origins of Sex Differences in Science‘. Mais surtout, les données sont téléchargeables ici.

> library(foreign)
> base=read.dta("http://www.ats.ucla.edu/stat/stata/examples/long/couart2.dta")
> head(base)
  art   fem     mar kid5  phd ment
1   0   Men Married    0 2.52    7
2   0 Women  Single    0 2.05    6
3   0 Women  Single    0 3.75    6
4   0   Men Married    1 1.18    3
5   0 Women  Single    0 3.75   26
6   0 Women Married    2 3.59    2

On dispose de plusieurs informations, recueillies auprès de 950 docteurs en biochimie. On a ainsi

  • le nombre de publications au cours des 3 dernières années de leur doctorat,
  • le sexe du docteur,
  • son statut marital (marié(e) ou célibataire), a priori à la fin du doctorat,
  • le nombre d’enfant(s) du docteur agé(s) de moins de 5 ans, là encore à la fin du doctorat
  • le nombre de publications du directeur de thèse (appelé ici mentor) au cours des 3 dernières années

(on a aussi le prestige du département, mais je trouve cette variable trop subjective). Comme on fait un comptage de publications, on pourrait espérer voir une loi de Poisson.

> n=nrow(base)
> lambda=mean(base$art)
> compte=table(base$art)
> N=as.numeric(names(compte))
> empirique=as.numeric(compte[as.character(0:19)])/n
> theorique=dpois(0:19,lambda)
> cbind(0:19,empirique,theorique)
        empirique theorique
 [1,]    0 0.3005  0.1839
 [2,]    1 0.2688    0.3114
 [3,]    2 0.1945    0.2636
 [4,]    3 0.0918    0.1487
 [5,]    4 0.0732  0.0629
 [6,]    5 0.0295  0.0213
 [7,]    6 0.0185  0.0060
 [8,]    7 0.0131  0.0014
 [9,]    8 0.0010  0.0003
[10,]    9 0.0021  0.0000
[11,]   10 0.0010  0.0000
[12,]   11 0.0010  0.0000
[13,]   12 0.0021  0.0000
[14,]   13     NA  0.0000
[15,]   14     NA  0.0000
[16,]   15     NA  0.0000
[17,]   16 0.0010  0.0000
[18,]   17     NA  0.0000
[19,]   18     NA  0.0000
[20,]   19 0.0010  0.0000
> plot(0:19,theorique,col="red",type="p")
> lines(0:19,empirique,type="p",col="blue ")

Mais non… graphiquement, on est relativement éloignés de la loi de Poisson

avec en rouge la fréquence que l’on aurait eu avec une loi de Poisson (de paramètre la moyenne obtenue sur l’ensemble de la population), et en bleu, la fréquence empirique.
Un test rapide permettra de convaincre les sceptiques,

> library(vcd)
> summary(goodfit(base$art,type="poisson"))

	 Goodness-of-fit test for poisson distribution

                      X^2 df     P(> X^2)
Likelihood Ratio 296.3715 13 1.381259e-55

Manifestement, il y a de l’hétérogénéité dans la population des doctorants. Il y a ceux qui publient beaucoup, ceux qui publient peu et ceux qui ne publient pas.

  • Le sexe et le nombre de publications

Peut-être qu’en discriminant par le sexe, on aurait, par sous-populations, des lois de Poisson. Si on étudie la distribution du nombre de publications pour les hommes

versus la distribution du nombre de publications pour les femmes

ça ne semble pas suffire pour expliquer toute l’hétérogénéité. Même s’il semble y avoir une différence ne serait-ce que sur le nombre moyen de publications par sexe,

> tapply(base$art,base$fem,mean)
     Men    Women 
1.882591 1.470309

Là encore, on peut faire un test rapide de comparaison de moyennes

> t.test(base$art~base$fem)

	Welch Two Sample t-test

data:  base$art by base$fem 
t = 3.3298, df = 885.945, p-value = 0.0009049
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.1692780 0.6552866 
sample estimates:
  mean in group Men mean in group Women 
           1.882591            1.470309

Visiblement les deux sous-populations sont différentes, mais le critère du sexe ne suffit pas à expliquer l’hétérogénéité.
Maintenant, notons que le critère retenu est le nombre de publications. La stratégie des personnes de sexe féminin peut être de viser moins de publications, mais dans de meilleures revues. En
tous les cas, le sexe pourrait expliquer en partie (et en partie seulement) l’hétérogénéité observée.

  • Le statut marital et le nombre de publications

On peut aussi comparer le statut marital, en particulier marié(e)

versus célibataire

Là encore, on peut comparer les moyennes, et faire un test d’égalité. Mais cette fois, l’influence est beaucoup plus faible que le sexe.

> t.test(base$art~base$mar)

	Welch Two Sample t-test

data:  base$art by base$mar 
t = -1.1869, df = 710.307, p-value = 0.2356
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.4033993  0.0994165 
sample estimates:
 mean in group Single mean in group Married 
             1.592233              1.744224

la différence entre les moyennes n’étant pas, ici, significative.

  • Les enfants et le nombre de publications

Une variable probablement plus intéressante que le statut marital pourrait être le nombre d’enfants (moi qui finissait la rédaction de ma thèse alors que ma deuxième apprenait à marcher On a la distribution suivante pour le nombre de publications pour un(e) docteur(e) ayant au moins un enfant,

versus un(e) docteur(e) sans enfant,

Mais la encore, le test de comparaison de moyenne n’est pas vraiment concluant,

> base$kid=base$kid5>0
> t.test(base$art~base$kid)

	Welch Two Sample t-test

data:  base$art by base$kid 
t = 0.6136, df = 646.359, p-value = 0.5397
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1803154  0.3442384 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.721202            1.639241

Certes, les docteur(e)s sans enfant publient plus, mais ce n’est pas non plus flagrant. On pourrait se dire que ce sont surtout les femmes qui ont eu des enfants qui sont pénalisées (la encore, car l’objectif est d’expliquer une quantité de publications), mais tout d’abord autant les hommes que les femmes répercutent le fait d’avoir eu, ou non, un enfant; mais surtout la variable de sexe est beaucoup plus discriminante que le fait d’avoir eu des enfants,

> baseH=base[base$fem=="Men",]
> baseF=base[base$fem=="Women",]
> t.test(baseH$art~baseH$kid)

	Welch Two Sample t-test

data:  baseH$art by baseH$kid 
t = 1.3831, df = 491.989, p-value = 0.1673
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1136025  0.6538074 
sample estimates:
mean in group FALSE  mean in group TRUE 
           2.011628            1.741525 

> t.test(baseF$art~baseF$kid)

	Welch Two Sample t-test

data:  baseF$art by baseF$kid 
t = 0.9525, df = 137.906, p-value = 0.3425
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1764049  0.5043375 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.501466            1.337500
  • Impact du directeur de thèse sur le nombre de publications

Une fois analysées les variables personnelles, si on regardait l’impact du directeur de thèse sur le travail de l’étudiant ? Qui sait, les variables pourraient être corrélées (je ne parle pas ici de relation causale). Par exemple, on peut regarder la distribution du nombre de publications pour un(e) docteur(e) ayant eu un directeur de thèse qui publiait peu (moins que la moyenne ici)

versus un directeur de thèse qui publiait beaucoup (plus que la moyenne),

Cette fois, on a une nette différence

> base$publis=base$ment>=9
> t.test(base$art~base$publis)

	Welch Two Sample t-test

data:  base$art by base$publis 
t = -6.7568, df = 497.017, p-value = 3.96e-11
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -1.2506429 -0.6871632 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.341338            2.310241

Un directeur de thèse qui publie semble avoir des étudiant(e)s qui publient aussi (ou réciproquement).

  • Régression de Poisson et nombre de publications

La population était initialement très hétérogène. Trop hétérogène pour une loi de Poisson, avec une variance valant plus du double de la moyenne empirique,

> var(base$art)/mean(base$art)
[1] 2.191358

> summary(glm(art~1,data=base,family=quasipoisson))

Call:
glm(formula = art ~ 1, family = quasipoisson, data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8401  -1.8401  -0.5770   0.2294   7.5677  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.52644    0.03761   14.00   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 2.191389)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1817.4  on 914  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Si on régresse sur les différentes variables on retrouve les interprétations que nous avions eues,

> summary(glm(art~mar+fem+I(kid5>0)+ment,data=base,family=quasipoisson))

Call:
glm(formula = art ~ mar + fem + I(kid5 > 0) + ment, family = quasipoisson, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.5080  -1.5615  -0.3626   0.5614   5.4494  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      0.338732   0.081482   4.157 3.53e-05 ***
marMarried       0.149527   0.085110   1.757  0.07928 .  
femWomen        -0.218643   0.074151  -2.949  0.00327 ** 
I(kid5 > 0)TRUE -0.250286   0.085834  -2.916  0.00363 ** 
ment             0.026000   0.002656   9.788  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 1.837376)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1641.1  on 910  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

On peut d’ailleurs faire une régression GAM avec une estimation par splines de l’effet du nombre de publications du directeur de thèse. Ce qui permet de légitimer le fait que la variable intervient de manière linéaire

On pourrait noter la concavité de la courbe qui semble sous-entendre qu’un directeur qui publie énormément peut toujours s’occuper de ses étudiants – et les faire publier plus qu’un directeur qui publie moins, mais en première approximation, la courbe peut être supposée linéaire.
On peut aussi regarder rapidement qui ne publie pas pendant sa thèse,

> summary(glm((art==0)~mar+fem+I(kid5>0)+ment,data=base,family=binomial))

Call:
glm(formula = (art == 0) ~ mar + fem + I(kid5 > 0) + ment, family = binomial, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1609  -0.9043  -0.7019   1.2697   2.5221  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -0.28843    0.18591  -1.551   0.1208    
marMarried      -0.33565    0.18800  -1.785   0.0742 .  
femWomen         0.23697    0.15829   1.497   0.1344    
I(kid5 > 0)TRUE  0.42957    0.19032   2.257   0.0240 *  
ment            -0.08141    0.01261  -6.455 1.08e-10 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1118.7  on 914  degrees of freedom
Residual deviance: 1052.0  on 910  degrees of freedom
AIC: 1062.0

Number of Fisher Scoring iterations: 5

Certain(e)s ont eu un enfant, et semblent avoir été pénalisé(e)s, mais si un doctorant n’a pas publié, c’est presque uniquement à cause de son directeur (qui publie peu).
Pour reprendre l’interprétation sur l’analyse du nombre de publications, par exemple, le fait d’être marié, ou pas, a peu d’influence sur le nombre de publications. Par contre les femmes – toutes choses étant égales par ailleurs – publient 20% de moins que les hommes,

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[3])
 femWomen 
0.8036087

Par contre le fait d’avoir – ou pas – des enfants semble ici beaucoup plus fort, car les docteur(e)s qui ont eu au moins un enfant publient 25% de moins que ceux qui n’en ont pas eu.

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[4])
I(kid5 > 0)TRUE 
      0.7785778

Aussi, une femme ayant eu un enfant aura publié 37% de moins qu’un homme qui n’a pas eu d’enfants. Toutes choses étant égales par ailleurs.

  • Toutes choses étant égales par ailleurs ?

C’est ici qu’est l’arnaque, car rien ne garantie que toutes choses puissent être égales par ailleurs. Par exemple, une femme avec enfant (ou pas) peut elle avoir un directeur de thèse qui publie beaucoup ? N’y a-t-il pas un peu d’endogénéité cachée ?
Si on regarde le nombre de publications du directeur de thèse versus le sexe de l’étudiant,

> t.test(base$ment~base$fem)

	Welch Two Sample t-test

data:  base$ment by base$fem 
t = 2.7101, df = 901.2, p-value = 0.006854
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.4587075 2.8673525 
sample estimates:
  mean in group Men mean in group Women 
           9.532389            7.869359

on note que les directeurs de thèse des hommes publient plus que les directeurs de thèse des femmes. Autrement dit, si les femmes publient moins, c’est peut être parce que les directeurs de thèses qui pourraient les aider à publier ne les prennent pas en thèse, non ? Car il n’y a pas de lien entre le fait d’avoir ou pas des enfants à la fin de sa thèse, et le nombre de publications de son directeur,

> t.test(base$ment~base$kid)

	Welch Two Sample t-test

data:  base$ment by base$kid 
t = -1.5626, df = 499.06, p-value = 0.1188
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.5485856  0.2905184 
sample estimates:
mean in group FALSE  mean in group TRUE 
           8.377295            9.506329

Bref, c’est manifestement le directeur de these qui influence le plus le nombre de publication d’un doctorant. Mais il serait intéressant d’avoir les memes statistiques 5 ans apres la soutenance…

Erratum to the paper on heatwave modelling

An erratum to the paper published on heatwave modelling should appear soon, on http://link.springer.com/article/…

In my article, On the return period of the 2003 heat wave, the sentence “summer of2003 was by far the hottest summer since 1500” is wrongly attributed to ‘Luterbacher et al. (2003)’. My probable intent was to cite Casty et al. (2005) which actually states “the years 1994, 2000, 2002 and particularly 2003 were the warmest since 1500”. Note that a similar statement can also be found in Stott et al. (2004), “The summer of 2003 was probably the hottest in Europe since at latest ad 1500”.

Regarding the statement (in reference to Luterbacher et al., which was actually published in 2004) that “their estimate of the return period of that event is 250 years”, the attribution again is incorrect. Such an estimate can be found in Brown et al. (2005), “best estimate of a 1 in 250 year event” or Lowe et al. (2004) “the same event is likely to occur much more frequently; most likely, once every 250 years”, among others.
I regret these inappropriate attributions and inadequate referencing.
  • Brown S, Stott P, Clark R, Temperature Extremes (2011) the Past and the Future. Hadley Centre for Climate Prediction and Research Working Paper. [online http://stabilisation.metoffice.com/…]
  • Casty C, Wanner H, Luterbacher J, Esper J, Böhm R (2005) Temperature and precipitation variability in the European Alps since 1500. Int J Climatol 25(14):1855–1880, 30 November 2005 CrossRef
  • Lowe J, Smith F, Jenkins G, Pope V (2004) Uncertainty, risk and dangerous climate change. Hadley Centre [online http://cedadocs.badc.rl.ac.uk/247/
  • Luterbacher J, Dietrich D, Xoplaki E, Grosjean M, Wanner H (2004) European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303:1499–1503 CrossRef
  • Stott PA, Stone DA, Allen MR (2004) Human contribution to the European heatwave of 2003. Nature 432:610–614 (2 December 2004) CrossRef

Loi des grands nombres et théorème central limite en assurance

Dans le dernier numéro de la revue Risques figure un petit article sur la loi des grands nombres et le théorème central limite en assurance (ici) Le but était de remettre un peu les pendules à l’heure sur la différence entre les implications de ces deux résultats, et surtout d’amorcer une discussion sur l’importance des hypothèses (en particulier celle de risques indépendants).

Beta kernel and transformed kernel

This Thursday I will give a talk at Laval University, on “Beta kernel and transformed kernel : applications to copula density estimation and quantile estimation“. This time, I will talk at the department of Mathematics and Statistics (13:30 at the pavillon Adrien-Pouliot). “Because copulas have bounded support (the unit square in dimension 2), standard kernel based estimators of densities are (multiplicatively) biased on borders and in corners of the support. Two techniques can be used to avoid that underestimation: Beta kernels and Transformed kernel. We will describe and discuss those two techniques in the first part of the talk. Then, we will see that it is possible to combine those two techniques to get nice estimator of several quantities (e.g. quantiles): transform the data to get on the unit interval – using a transformed kernel – then estimate the (transformed) quantile on [0,1] using a beta kernel, then get back on the initial support. As we will see on simulations, that technique can be better than standard quantile estimators, especially when data are heavy tailed.” Slides can be downloaded here.

  • kernel based density estimation

Kernel based estimation are a popular (and natural) technique to estimate densities.  It is simply and extension of the moving histogram:

so we count how many observations are a the neighborhood of the point where we want to estimate the density of the distribution. Then it is natural so consider a smoothing function, i.e. instead of a step function (either observations are close enough, or not), it is possible to give weights to observations, which will be a decreasing function of the distance,

With a smooth kernel, we have a smooth estimation of the density

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-01.gif

Then it is possible to play on the bandwidth, either to get a more accurate estimation of the density, but not that smooth (small bias but large variance),

or a smoother one (large bias, but small variance),

In R, it is simply

> X=rnorm(100)
> (D=density(X))
 
Call:
	density.default(x = X)
 
Data: X (100 obs.);	Bandwidth 'bw' = 0.3548
 
       x                   y            
 Min.   :-3.910799   Min.   :0.0001265  
 1st Qu.:-1.959098   1st Qu.:0.0108900  
 Median :-0.007397   Median :0.0513358  
 Mean   :-0.007397   Mean   :0.1279645  
 3rd Qu.: 1.944303   3rd Qu.:0.2641952  
 Max.   : 3.896004   Max.   :0.3828215  
 
> plot(D$x,D$y)
  • Beta kernel

The idea of Beta kernel is to consider kernels having support [0,1]. In the univariate case,

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-06.gif

where http://freakonometrics.blog.free.fr/public/perso3/kernel-f-07.gif is the density of a Beta distribution, i.e.

http://freakonometrics.blog.free.fr<br />
/public/perso3/beta-distribution.gif

For additional material, I have uploaded some R code to fit copula densities using beta kernels,

library(copula)
beta.kernel.copula.surface = function (u,v,bx,by,p) {
s = seq(1/p, len=(p-1), by=1/p)
mat = matrix(0,nrow = p-1, ncol = p-1)
for (i in 1:(p-1)) {
a = s[i]
for (j in 1:(p-1)) {
b = s[j]
mat[i,j] = sum(dbeta(a,u/bx,(1-u)/bx) *
dbeta(b,v/by,(1-v)/by)) / length(u)
} }
return(data.matrix(mat)) }

Then we can used it to see what we get on a simulated sample

library(copula)
COPULA = frankCopula(param=5, dim = 2)
X = rcopula(n=1000,COPULA)
p0 = 26
Z= beta.kernel.copula.surface(X[,1],X[,2],bx=.01,by=.01,p=p0)
u = seq(1/p0, len=(p0-1), by=1/p0)
persp(u,u,Z,theta=30,col="green",shade=TRUE,
box=FALSE,zlim=c(0,6))

http://freakonometrics.free.fr/copula-kernel-beta.gif
(yes, the surface is changing… to illustrate the impact of the bandwidth on the estimation).

  • transformed kernel estimation

I the talk, I will also mention the transformed Kernel estimate, as introduced in the book on L1 density estimation by Luc Devroye and Laszlo Györfi (the book can be downloaded here). I probably spend a few minutes on the original chapter, in order to provide another application of that techniques (not only to estimate copula densities, but here to estimate quantiles of heavy tailed distribution). In the univariate case, the R code is the following (here I consider two transformation, the quantile function of the Gaussian distribution, and the quantile function of the Student distribution with 3 degrees of freedom),

set.seed(1)
sample=rbeta(100,4,3)
 
transfN = function(x){
Y=qnorm(sample)
f=density(Y,from=-4,to=4,n=2001)
ny=sum(f$x<=qnorm(x)); 
  g=f$y[ny]/dnorm(qnorm(x))
return(g)
}
 
df0=3
 
transfT = function(x){
Y=qt(sample,df=df0)
f=density(Y,from=-4,to=4,n=2001)
ny=sum(f$x<=qt(x,3)); 
  g=f$y[ny]/dt(qt(x,df=df0),df=df0)
return(g)
}
 
tN=Vectorize(transfN)
tT=Vectorize(transfT)
 
u=seq(.01,.99,by=.01)
vN=tN(u)
vT=tT(u)
plot(u,vN,type="l",lwd=3,col="blue")
lines(u,vT,lwd=3,col="green")
lines(u,dbeta(u,4,3),col="red",lty=2)

The density estimation is the following,

(the red dotted line is the true density, since we work on a simulated sample). Now, let us get back on the initial chapter,

In the book, this is introduced as follows,

The original idea we add it to use this kernel based estimator for copulas, i.e. since we can estimate densities in high dimension with unbounded support, using

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-02.gif

the idea is to transform marginal observations,

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-10.gif

and to use the fact that the associated copula density can be written

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-12.gif

to derive an intuitive estimator for the copula density

http://freakonometrics.blog.free.fr/public/perso3/kernel-f-13.gif

An important issue is how do we choose the transformation

And Luc Devroye and Laszlo Györfi mention that this can be used to deal with extremes.

well, extremes are introduced through bumps (which is not the way I would have been dealing with extremes)

and note that several results can be derived on those bumps,

e.g.

Then, there is an interesting discussion about estimating the optimal transformation

and I will prove that this can be an extremely interesting idea, for instance to estimate quantiles of heavy tailed distribution, if we use also the beta kernel estimator on the unit interval. This idea was developed in a paper with Abder Oulidi, online here.

Remark: actually, in the book, an additional reference is mentioned,

but I have never been able to find a copy… if anyone has one, I’d be glad to read it…

Talk at Laval University at the Actuarial Seminar

I was last Friday at Laval University for a conference by David Cummins and Mary Weiss (here). I will be back tomorrow, this time to give a talk, on “distorting probabilities in actuarial science” (the talk will be extremely close to the one I gave at McGill in November). “In this talk, we will first get back on properties of distortion operators for pricing financial and insurance risks. Based on the dual version of the expected utility framework, we will see how distorted risk measures have been introduced, from VaR and TVaR, to Esscher premium and Wang’s measures. Then we will discuss extensions in higher dimension. We will discuss tail properties of distorted copulas (in the particular case of Archimedean copulas). A natural application will be aging problems (in survival analysis or in credit risk).” Slides can be downloaded from here.

 

This talk can be seen as a first part, the second one behing the talk I will give in 15 days, again at Laval University, but this time for the Seminar of Statistics. The talk will be on “Beta kernel and transformed kernel : applications to quantile estimation, and copula density estimation“.

Talk at Desjardins General Insurance

This afternoon, I will give a talk at the seminar of the R&D department at Desjardins General Insurance, on correlation in claims reserving. A lot of interesting papers have been published recently on that topic. On multivariate Chain Ladder, some interesting articles have been published, e.g. the one by Carsten Prohl and Klaus Schmidt (here) or the one by Michael Merz and Mario Wuthrich (there).

But I think another interesting perspective (so far, not in claims reserving, but one should find some time to look at it) should be about multivariate regression (multivariate GLM’s), e.g.

All that will be mentioned in the talk. Slides can be downloaded here,

The dataset used in the example can be obtained with the code below

> P.corp=read.table("http://freakonometrics.blog.free.fr/public/data/auto-corporel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.corp=as.matrix(P.corp) > n=nrow(P.corp) > P.mat =read.table("http://freakonometrics.blog.free.fr/public/data/auto-materiel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.mat=as.matrix(P.mat) > P.mat=P.mat[1:n,1:n] >  P.mat = P.mat[2:10,1:9] >  P.corp= P.corp[2:10,1:9] > n=9 > P.tot = P.mat + P.corp

When will my papers appear as references (if they do…) ?

Following my post on citations in academic journals, I wanted to go one step further in the understanding of the dynamic of citations. So here, the dataset looks like that: for each article, we have the name of the journal, the year of publication (also the title of the article, but here we do not use it, as well as the authors), and more interesting, the number of citations in journals (any kind of academic journal) published in 1996, 1997, …, 2011. Of course, articles published in 1999 might have their first citation only starting in 1999.

base[1000:1002,]
     Publication.Year
7188             1999
7191             1999
7195             1999
     Document.Title
7188 Sequential inspection 
7191 On equitable resource approach
7195 Method for strategic  
                                        Authors     ISSN       Journal.Title
7188                         Yao D.D., Zheng S. 0030364X Operations Research
7191                                    Luss H. 0030364X Operations Research
7195 Seshadri S., Khanna A., Harche F., Wyle R. 0030364X Operations Research
     Volume Issue X139 DEV1996 DEV1997 DEV1998 DEV1999 DEV2000 DEV2001 DEV2002
7188     47     3    0       0       0       0       0       1       0       2
7191     47     3    0       0       0       0       0       0       2       0
7195     47     3    0       0       0       0       0       0       0       0
     DEV2003 DEV2004 DEV2005 DEV2006 DEV2007 DEV2008 DEV2009 DEV2010 DEV2011
7188       0       0       0       1       0       0       0       0       0
7191       3       4       1       4       4       8       4       6       1
7195       0       1       2       2       1       0       1       0       0
     X130655 X0 X130794
7188       4  0       4
7191      37  0      37
7195       7  0       7

The first step is to aggregate data, not to look at each article, but to look at all paper published in 1999 (say). And then, we look at the number in citations the year of publication, the year after, two years after, etc. It will appear in a triangle since if we look at articles published in 2010, there is only on possible year for citations (2010, since I removed 2011).

VOL=rev(unique(base$Volume))
VOL=VOL[is.na(VOL)==FALSE]
TRIANGLE=matrix(NA,16,16)
for(v in VOL){
k=k+1
sb=base[base$Volume==v,9:24]
sb=sb[is.na(sb[,1])==FALSE,]
TRIANGLE[k,1:(17-k)]=apply(sb,2,sum)[k:16]}

Then, a standard idea (at least in insurance business, for claims payment development) is to consider that data are Poisson distributed, and the number of citations should depend on the year of publication of the article (a row effect) and the development (how many years after are we looking at, i.e. a column effect). More formally, let http://freakonometrics.blog.free.fr/public/perso2/citationD01.gif denote the number of citations of articles published year http://freakonometrics.blog.free.fr/public/perso2/citD02.gifduring year http://freakonometrics.blog.free.fr/public/perso2/citD03.gif (or after http://freakonometrics.blog.free.fr/public/perso2/citD04.gif years). And we assume that http://freakonometrics.blog.free.fr/public/perso2/citD05.gif

TRIANGLE=TRIANGLE[-16,]
TRIANGLE=TRIANGLE[,-16]
Y=as.vector(TRIANGLE)
YEAR=rep(1996:2010,15)
DEV =rep(1:15,each=15)
baseT=data.frame(Y,YEAR,DEV)
reg=glm(Y~as.factor(YEAR)+as.factor(DEV),
data=baseT,family=poisson)

Since those are incremental values, in order to look at the paper of distribution, we need to sum them on a line. Thus, we can plot

http://freakonometrics.blog.free.fr/public/maths/dev-cl-biblio-1.gif
http://freakonometrics.blog.free.fr/public/maths/dev-cl-biblio-2.gif

(because we used factors, the first component has been replaced by the constant in the regression) or a normalized version to compare among journals. For instance, we would like to get 100 citations over 15 years.

DYN=exp(c(reg$coefficients[1],reg$coefficients[1]+
reg$coefficients[16:29]))
DYNN=cumsum(DYN)/sum(DYN)
plot(0:15,DYNN)

And this is what we get, for several academic journals,

The pattern is rather different. For instance, in Health Economics, citations is a quick process: more than 40% of citations obtained over 15 years, were obtained during the first 4 years. On the other hand, in the Journal of Finance, it is much smaller: less than 15% of the citations were obtained during the first 4 years (on average). So it means that comparing citation based index (namely g or h) is a difficult exercise, especially with you researchers in different areas. The same gor index for young researcher, publishing either in Stochastic Processes and their Applications or Annals of Statistics, means that after 3 years, it can be 50% higher.


Now it is possible to look more into details, with below JRSS-B (on applied statistics). Note that here, citations come extremely slowly… to it might not be a good “strategy” (assuming that a researcher’s target is simply to get – quickly – a high citation index) for a young researcher to publish in JRSS-B

On the other hand, Biometrika is much faster (both are on applied statistics, but we’ve seen here that they were not in the same cluster)

We can also observe that Annals of Probability
and Stochastic Processes and their Applications

have (almost) similar patterns (SPA might be a bit faster). Anyway, I have been surprised to see that in theoretical journals citations are extremely fast. Especially if we compare with the Journal of Finance for instance

where I though citations were extremely fast. But I might have a non-correct interpretation: it might simply mean that in the Journal of Finance it is common to cite old papers (published 10 or 15 years ago), maybe more common that in stochastic processes…
Anyway, all suggestions about the interpretation are welcomed !

Talk at Laval University on natural catastrophes

On Tuesday, I will be giving a talk at the Département de finance, assurance et immobilier, at the Faculté des sciences de l’administration. The talk will be on natural catastrophes, and on government intervention. The slides will be upladed soon (since we are still revising the paper we wrote Benoît, cf here: actually, we did not look at EU maximizers but RDEU maximizer with a quantile-based distortion). I will write a more detailed post once the working paper is finished.

Going further on journal clustering: looking within a cluster

Following my post on academic journals, Miss Lambda asked me about journals in a very specific area, namely agricultural, environmental and energy economics. I found it interesting since I do not have any idea about journal that can be in that domain of research. So I looked at some journals form the French CNRS list (online here). Hence, I have been looking for words in the title of 26,000 articles, published in 29 journals, in Agricultural, Environmental and Energy Economics. I considered American Journal of Agricultural Economics (AJEE), Ecological Economics (EE), Journal of Environmental and Economic Management (JEEM), Climate Policy (CP), Energy Economics (EE), Energy Journal (EJ), Energy Policy (EP), Environment and Planning A, B, C, D (EP-A, B, C, D), Environmental and Resource Economics (ERE), Environmental Modeling and Assessment (EMA), European Review of Agricultural Economics (ERAE), Resource and Energy Economics (REE), Agricultural Economics (AE), AMBIO: A Journal of the Human Environment (AMBIO), Australian Journal of Agricultural and Resource Economics (AJARE), Canadian Journal of Agricultural Economics (CJAE), Climatic Change (CC), Ecological Modeling (EM), Energy Studies Review (ESR), Environment and Development Economics (EDE), Environmental Science and Policy (ESP), Environmental Values (EV), Food Policy (FP),
Global Environmental Change (GEC), Journal of Agricultural Economics (JAE), Society and Natural Resources (SNR), Water Resources Research (WR).

Now if we look at the principal component analysis, and projection of the journals on the first two axis, we have

On that graph, it is hard to say anything… The only thing is see is that on the lower part, we have journals focusing on modeling issues.

The top of the most common words in those journals is the following,

> colnames(MATRICE[,1:32])
[1] "climate"       "analysis"      "environmental" "change"
[5] "energy"        "model"         "policy"        "case"
[9] "water"         "economic"      "management"    "food"
[13] "study"         "development"   "market"        "agricultural"
[17] "approach"      "effects"       "carbon"        "global"
[21] "impact"        "production"    "land"          "china"
[25] "assessment"    "forest"        "impacts"       "urban"
[29] "demand"        "spatial"       "emissions"     "modeling"

If we look at their projections on the first two axis, we have (it looks like the second axis has been inverted here)

i=or if we focus on the top 30

Now, if we look at clusters, and use a hierarchical model, we obtain

So here, a dozen journals are extremely close. They are more focusing on agricultural issues. Note that Energy Economics and Energy Policy are in the same cluster, but quite far away from Energy Journal (which looks strange). We can also observe that Climatic Change alone, far away from all the other journals (actually, it is a journal were I just got a paper accepted… and since my areas of research are quite far away from agricultural economics, I can understand that).

Think academic journals look the same ? Well, some do…

We have seen yesterday that finding an optimal strategy to publish is not that simple. And actually, it can be even more difficult in the case the journal rejects the paper (not because it is not correct, but because “it does not fit” with the standards, the quality of the journal, the audience, the editor’s mood, or whatever). The author has basically two choices,

  • forget about the article and move to something else (e.g. start a blog where he/she will be the author and the editor)
  • pretend that the article is worth publishing and then try to find another journal with similar interests


But this last choice is not that easy, since sometimes the author think that this journal was indeed the one that should publish it (e.g. all the articles on the subject have been published in that journal).
So I was wondering if there were clusters of journals, i.e. journals that publish almost the same kind of articles (so that next time one of my paper is rejected by the editor, I just go to for some journal in the same cluster).
So what I did is extremely simple: I looked at articles titles and looked for correlations between words frequency (I could have done that in key words, but I am not a big fan of those key words). I looked at 35 journals (that are somehow related to my areas of interest) and looked at titles of all articles published over the last 20 years. Then I kept the top 1000 of words, and I removed standard short words (“a“, “the“, “is“, etc). Actually, my top words looks like

"models" "model" "data" "estimation" "analysis" "time" 
"processes" "risk" "random" "stochastic" "regression" 
"market" "approach" "optimal" "based" "information" 
"evidence" "linear" "games" "bayesian" "theory" "effects"
"distribution" "multivariate" "tests" "markets" "markov"
"equilibrium" "dynamic" "process" "distributions" 
"application" "stock" "likelihood"

Then, I ran a principal component analysis on my dataset (containing 960 variables – here words – and 35 observations – here journal names).

library("FactoMineR")
res.pca = PCA(MATRICE, scale.unit=TRUE, ncp=5, 
graph=FALSE)
plot.PCA(res.pca, axes=c(1, 2), choix="ind")

The projection of the journals on the first two axis looks like that

Here, we can clearly observe some clusters : on the up-left Journal of Finance and Journal of Banking and Finance (say financial journals) on the top-right Biometrika, Biometrics, Computational Statistics and Data Analysis and Journal of Econometrics (JASA is not far away, i.e. applied statistics journal). And below, on the right, Stochastic Processes and their Applications, Annals of Applied Probability, Journal of Applied Probability, Annals of Probability, Proceedings of AMS and Topology and Applications (ie more theoretical journal).
Note that the projection is rather robust: if I consider my first 200 words, the graph is the same

In order to go further in the interpretation, we can also plot variables, i.e. words from titles,

where we cannot distinguish anything. So if I just look at my top 30, here they are,

On top left we see market(s), risk or information; on top right analysis, effects, models or tests; while below we see Markov or process(es). And we can observe interesting facts: in finance in statistics, we talk about dynamics while in theoretical (mathematical) journal it is about processes.
But the goal was to find cluster, i.e. classes of journals that publish papers with similar titles.

DISTANCE = dist(MATRICE)
cah = hclust(DISTANCE) 
plot(cah)

Here we have

If some classes a rather natural (Journal of Applied Proba. and Advances in Applied Proba.or Economic Theory, Journal of Economic Theory and Journal of Mathematical Economics) some strong correlation are not simple to understand, (e.g. Insurance: Mathematics and Economics and Management Science or Annals of Statistics and the Journal of Multivariate Analysis).
Again, it might be possible to spend hours on the graphs, but if I want – someday – to submit something to one of those journals, I guess I have to stop here, and move to something else…

 

The longer the better ?

No, this time it is not a post for Valentine’s day… It is simply that a few days ago, on Gaïa Universitas, Rachel mentioned a study on academic journals in astronomy. The study is online here, written by Krzysztof Zbigniew Stanek. « Naively, one would expect longer papers to have larger impact (i.e., to be cited more) – how long a paper should be to maximize its impact? Is it better to write several shorter papers or one longer paper? ». Actually, I did not expect that, and I was truly surprised to see that the number of citations was an increasing function of the number of pages. So I had a look at several academic in economics and mathematics. The first journal I looked at is the Journal of Finance. Here, the relationship between the number of pages and the number of citation is strong: a 40 page paper is – on average – two times more cited than a 20 page paper. Actually, the shape of the regression curve is rather close to the one obtained in Krzysztof’s paper,

The colors are due to the fact that some articles were published in 2010, and some in 1995. Obviously, the number of citations is a function of the year of publication (a longer post will be published soon on the dynamics of the citation process). I considered only articles published before 2005 to fit a regression model (here a spline regression, to see if the function is linear). And indeed, a longer paper has more chance to be cited.
Then, I looked at the Journal of Multivariate Analysis, which is a mathematical journal, where (theoretical) econometricians can publish.

Here the pattern is rather different: if we do not take into account short notes, the number of citations is independent of the size of the paper.
The more I look at those graphs, the more disturbing I find them…

  • on the one hand, I strongly believe that the size of the paper has nothing to do with the quality (or the importance) of the paper: so in some sense, I was expecting a flat regression, like the one we see above, with the Journal of Multivariate Analysis. It might come from the fact that in statistics (for instance) with empirical processes, e.g., proofs are extremely long and technical (and papers are long), while on stochastic orderings, proofs are extremely short (and papers rather short). But both can appear in the same journal…
  • on the other hand, researchers are more and more evaluated, and a common tool is to look at citation indexes (the called « publish or perish » paradigm). The more citations, the better the researcher, something like that… So sometimes, when we have a quite long paper, the question that naturally arises is: why not splitting the paper in two ? with flat regression, it is, indeed, optimal to split the paper (if possible) in two, since two papers with 20 page each might yield two times more citations than one. But with a curve like the one we see with the Journal of Finance, such a split is not an issue…

Anyway, I did really like the conclusion of Krzysztof’s paper « This paper will not be submitted to any journal, but please feel free to cite it as often as possible, or better yet cite my regular astronomical papers »….

Uncertainty in claims reserving (WIM)

I will be talking on Friday, at the 1st Québec-Ontario Workshop on Insurance Mathematics (so called WIN, already mentioned here). The program is now online here. My talk will be on Solvency II requirements in claims reserving (slides can be found here). Even if Canada will not adapt the Solvency II European capital test, it looks like Solvency II will matter to Canadian insurance companies (as mentioned here).

The return period of the 2003 heat wave

The paper on “on the return period of the 2003 heat wave” has been published online here, before Christmas. It should be published soon by Climatic Change, http://link.springer.com/article/...

Extremal events are difficult to model since it is difficult to characterize formally those events. The 2003 heat wave in Europe was not characterized by very high temperatures, but mainly the fact that night temperature were no cool enough for a long period of time. Hence, simulation of several models (either with heavy tailed noise or long range dependence) yield different estimations for the return period of that extremal event.

To go further on the impact on mortality (which was not the aim of the paper), there is a paper in Nature (here).