# Exchange rates and correlation matrices

I wanted to upload here a small problem I started to work on…. unfortunatley, I could not find (yet) a proper answer. Any comments and suggestions are welcomed. The problem is simple: consider 3 random variables that can be interpreted as exchange rates. Thus X is the exchange rate USD versus EUR, Y is the exchange rate EUR versus YEN, while Z is the exchange rate USD versus YEN. We then have the following constraint XY=Z. Given that constraint, what are the implied constraints on the correlation matrix of random vector (X,Y,Z) ?

• Bounds on correlation

There are classically standard constraints on correlation. For instance, from Cauchy-Schwarz’s inequality

but those bounds are not necessarily sharp, given distributions of  and . From Fréchet-Hoeffding bounds, we know that

where . Further, correlation matrices are necessarily positive-semidefinite matrices.
To go further to find possible bounds, recall that

and therefore

i.e.

Similarity,  will be a function of .
Thus, given marginal distributions of  and , covariances (and also correlations  and ) are functions of , and .
Somehow, we have obtained 2 degrees of freedom: if there are no additional constraints on those two coefficients, there should be no other constraints than having positive-semidefinite matrices.
I guess that  and  should not be connected. Those quantities are related to the co-skewness* and they should be different when dependence is asymmetric.
Given a 3 parameter copula for pair , respectively describing global dependence (that can be related to  and tail asymmetry (that should be related either to  and depending on the corner) it should be possible to derive any kind of triplets.

• Running simulations

I have tried to generate a lot of samples. The idea is to draw pairs from mixtures of Marshall-Olkins’ copula (wich are asymmetric), where parameters where generated randomly. I draw 250 values according to that copula. Then, the third variable is obtained taking the product of normalized pairs (or Studentized). And then calculate the associated correlation matrix.

1.  draw  from  (but actually any positive distribution should work) and
2. do 250 times the following (index )
• draw  with a Bernoulli  distribution
• draw&nbs
p; from copula  if  and  if
• set  and , and
1. estimate the correlation of sample

In dimension 3, I have the following graph,

In order to get a better understanding, I look at slices, for instance when . On the graph on the right, to dark area are values that can not be obtained since correlation have to be positive-semidefinite. Dots are points that have been obtained in a particular scenario. I have also included the convex hull if the points obtained with all the scenarios (here 50,000 scenarios), since I assume that the set of positive correlation is necessarily convex. Graphs below have been obtained when  and

So far I have not been able to obtain all possible correlations, I guess…  If anyone has suggestions…
* Recall that the co-skewness matrix, in dimension 2, for a pair  is simply

where .

# Les tables de mortalité

Chose promise, chose dûe, un court billet expliquant les principales tables utilisées,

• Les tables TD et TV 88-90

Ces tables datent un peu, et si je continue à en parler en cours, c’est parce qu’elles sont simples à utiliser (et pour continuer à me faire croire que je n’ai pas vieilli depuis mes études). Cette table est d’ailleurs tellement sérieuse qu’on la retrouve dans la loi (ici), dans un arrêté d’avril 1993. La table dite TD 88-90 (pour Décès) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population d’hommes. Elle était appliquée pour le calcul des primes des contrats d’assurance décès. La table dite TV 88-90 (pour Vie) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population de femmes. Elle était appliquée pour le calcul des primes des contrats d’assurance en cas de vie. Ces tables peuvent se récupérer à l’aide des codes suivants,

```> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",

Ces tables ont été remplacées par les tables dites TH et TF, respectivement.

• Les tables TH et TF 00-02

Ces tables ont été établies à partir des données INSEE de la population française entre 2000 et 2002 et ont été lissées. Ce sont des tables générationnelles, qui nécessitent un correctif d’âge pour tenir compte des écarts de mortalité entre générations. Elles sont applicables à partir du 1er janvier 2006. L’institut des actuaires a proposé une “notice d’utilisation” en ligne ici, et Cimon en avait parlé sur son blog ().

• Les tables TPRV 95

La table TPRV 93 (pour Table Prospective de Rente Viagère) est un extrait de la table dite plancher pour la tarification des contrats de rente viagère. Elle a été publiée par l’arrêté du 28 juillet 1993 (ici sans les annexes), et correspond à une table prospective qui retrace la mortalité des générations 1887 à 1993 (les tables prospectives sont au programme de Master 2).

La TPRV 93 représente la table complète de la génération 1950. La table est en ligne ici (en csv) lisible sous R avec le code suivant,

```> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TPRV.csv",

Mais ce n’est que la première étape. On utilise ensuite un décallage d’âge.

On a parlé de ce point en TD (mais sur les taux de mortalité), c’est aussi ce qui s’appelle l’hypothèse de Rueff, et traduit la diminution des taux de mortalité par un rajeunissement (i.e. un gain sur les âges), i.e.

Notons qu’en fait,  dépend du niveau des taux d’intérêt. L’ensemble des tables de décallage sont en ligne ici.

• Les tables TGH et TGF 05

Il s’agit ici, là encore, d’une population de rentiers, comme pour la (ici pour les 38 pages du code). Aussi, ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici.

# from two to three…

A short post to give more details about the final remark in the course of Financial Econometrics, and more precisely the formula that can be found in the book of Philip Jorion,

Note that this formula can be found (perhaps written with slight changes) in several papers, e.g. the following sentence (on the http://www.bis.org/website),

or the following formula, on documents from the Bank of England website,

I recently pulished (in French, here) a paper on the Value-at-Risk, including the following graph,

Usually, three times the average over 60 trading days is the larger component, but during the financial crisis, it turned out that the daily component was almost three times higher than the average value over the past the months (this fact was mention by Paul Embrechts in some conference in Paris on risk measures).
The interpreation of the multiplicative k coefficient (which is from 2 to 3 in some publications, or which exceeds 3 in others) has been proposed in a paper of Gerhard Stahl, entitled three cheers. The idea is to use the Bienaymé-Tchebychev inequality. For random variables with finite variance, then

Recall that this inequality is simply a corrolary of Markov’s inequality

or for any increasing function

(taking function , applied to ). This upper bound can be far away from the true probability, see e.g. the gaussian case below, i.e. if  ,

```> z = seq(0,3,by=.01)
> P = 2*dnorm(k)
> U = 1/z^2
> plot(z,P,type="l",lwd=2,col="red",xlab="",ylab="")```

The ratio between the two is given below,

`> plot(z,U/P,type="l",lwd=2,col="purple",xlab="",ylab="",ylim=c(0,10))`

Note that it is possible to interprete the axis values as probabilities values, taking quantiles of the gaussian distribution

```> plot(pnorm(z),U/P,type="l",lwd=2,col="purple",xlab="",
+ ylab="",ylim=c(0,10),xlim=c(.9,1))
> abline(h=3,lty=2)```

The interpretation is that the upper bound is 3 times higher than the true probability in the Gaussian case when z is the quantile of the  distribution associated with probability level 99%.
Note that

• if z is the 95% quantile of the mathcal{N}(0,1) distribution, the ratio is 2 (1.92)
• if z is the 99% quantile of the mathcal{N}(0,1) distribution, the ratio is 3 (3.04)
• if z is the 99.55% quantile of the mathcal{N}(0,1) distribution, the ratio is almost 4 (3.88)
• if z is the 99.75% quantile of the mathcal{N}(0,1) distribution, the ratio is 5 (5.04)

A more formal explaination is to assume that X is symmetric, and then

Thus, if , i.e. , we have an upper bound for the  Value-at-Risk,

where the upper bound is the upper bound for the  Value-at-Risk for any distribution with finite variance and centred.
If  , then , i.e. .  But since,  for a  distribution, then

and further

# Les modèles en réassurance

Publication d’un papier sur les modèles en réassurance dans la revue Risques (ici), suite à un questionnement sur la pertinence des modèles classiques utilisés par les réassureurs. La question initiale était partie de la constatation – que l’on retrouve ici ou – sur l’utilisation (ou la mauvaise utilisation) de modèles sophistiques en finance de marché, en essayant d’expliquer que – d’un point de vue épistémologique au moins – les modèles des réassureurs étaient plus robustes. En particulier, on notera que les modèles les plus anciens utilisés par les réassureurs (en particulier la loi de Pareto comme je l’avais évoqué ici) ont eu une légitimité pratique durant plusieurs décennies avant d’être justifiés par la théorie des valeurs extrêmes. Je ferais d’ailleurs bientôt un billet sur l’histoire des valeurs extrêmes en statistiques, en revenant en particulier sur les travaux de Gumbel ou de Fréchet.

Sinon le code utilisé dans le papier est en ligne ici, et la base

```> sinpe = read.table("https://perso.univ-rennes1.fr/arthur.charpentier/sinpe.csv",header=TRUE,sep=";")
DSUR  MNTPE
1 19850206 240439
2 19851228 125674
3 19850504 488331
4 19851118 457347
5 19850220 990919
6 19851214 182939
> annee=as.numeric(substr(as.character(sinpe\$DSUR),1,4))
> sinistres=sinpe\$MNTPE[annee>1992]
> XS=sinistres/100000```

On se limite ici aux sinistres survenus après 1992. Si l’on visualise ces montants de sinistres, on obtient

```> datesur=as.Date(as.character(sinpe\$DSUR),"%Y%m%d")
> jour=datesur[annee>1992]
> plot(jour,sinistres/100000,xlab="",ylab="Coût individuel",cex=.5,ylim=c(0,600))
> ded=50
> abline(h=ded)```

On peut aussi faire le graphique Pareto-log-log,

```> library(evir)
> n=length(X)
> plot(log(sort(X)),log((n:1)/(n+1)),
+ xlab="Coûts des sinistres (logarithme)",ylab="Fonction de survie (logarithme)",cex=.8)
> out <- gpd(XS, 15)
> XI=as.numeric(out\$par.ests[1]); BETA=as.numeric(out\$par.ests[2])
> x0=seq(2,8,.01)
> lines(x,-1/XI*(x-log(15)),col="red")```

(l’ajustement de la loi de Pareto permettant de tracer la droite rouge) ou encore visualiser l’estimateur de Hill de l’indice de queue,

`> hill(X)`

Pour finir, le code suivant permet de calculer la prime pure (ou plus généralement une prime de Wang) soit de manière non-paramétrique (burning cost) ou bien en utilisant l’ajustement d’une loi de Pareto. Le papier étant un papier de vulgarisation, j’ai pris le seuil de manière arbitraire, sans aucune recherche de valeur “optimale”

```> DEDUC = seq(10,50,by=5)
> lambda=0;  seuil=5
> WG1=WG2=rep(NA,length(DEDUC))
> for(k in 1:length(DEDUC)){
+ deductible=DEDUC[k]
+ out <- gpd(XS, seuil)
+ XI=as.numeric(out\$par.ests[1]); BETA=as.numeric(out\$par.ests[2])
+ G0=function(x){1-pgpd(x+seuil, xi = XI, mu = seuil, beta = BETA)}
+ G=function(x){(G0(x+deductible-seuil))/(G0(deductible-seuil))}
+ F=function(x){pnorm(qnorm(G(x))+lambda)}
+ (wang1=integrate(F, 0, Inf))
+ X=XS[XS>deductible]
+ n=length(X)
+ FS= function(z){
+ m=rep(NA,length(z))
+ for(i in 1:length(m)){
+ m[i]=sum(X>z[i]+deductible)/n}
+ return(m)
+ }
+ G=function(x){pnorm(qnorm(FS(x))+lambda)}
+ (wang2=sum(G(seq(0,800,.01))*.01))
+ WG1[k]=as.numeric(wang1\$value)
+ WG2[k]=wang2
+ }
> plot(DEDUC,WG2-DEDUC,type='b',xlab="Niveau de la priorité ('00 000 euros)",ylab="Prime pure par sinistres réassuré",ylim=c(0,50))
> lines(DEDUC,WG1-DEDUC,type='b',col="red",pch=4)
> legend(10,50,c("Aujustement d'une loi de Pareto","'Burning cost'"),
+ col=c("red","black"),lty=1,cex=.8,pch=c(4,1))```

# Nikkei’s past experience vs. SP500 (in euros)

Following Michael’s idea (here), I wanted to go further, based on his intuition (and dataset that he kindly sent me, there). If we consider the two series of Nikkei index and SP500 index in euros, we have to following graph,

the code is simply the following (the merging function is simply here to avoid problem with different trading days: since we look at the index and not the return, it is the simplest way to deal with it).

```> library(RODBC)
> base = odbcConnectExcel(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/spx_nky_eurusd.xls",
> series1 = sqlQuery(base,query="select * from [Tabelle1\$A2:B8837]") # SPX
> series2 = sqlQuery(base,query="select * from [Tabelle1\$D2:E8631]") # NKY
> series3 = sqlQuery(base,query="select * from [Tabelle1\$G2:H8945]") # EURUSD
> odbcCloseAll()
> series4=merge(series1,series3)
> series4\$SPEUR=series4\$SPX/series4\$EURUSD
> series5=merge(series4,series2)
> x=(as.Date(series5[,1])-as.Date("01/01/0000","%d/%m/%Y"))/365.25
> yl=range(series5[,4])
> xl=c(1975,2010)
> plot(x,series5[,4],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="red",xlim=xl,ylim=yl)
> axis(1)
> axis(2, col="red")
> par(new=TRUE)
> yl=range(series5[,5])
> plot(x,series5[,5],axes=FALSE,xlab="",ylab="",type="l",
+ lwd=3,col="blue",xlim=xl,ylim=yl)
> axis(4, col="blue")
> mtext("SP500 in Euros", 2, line=2, col="red", cex=1.2)
> mtext("NKY", 4, line=2, col="blue", cex=1.2)```

Those two series series seem to have a similar pattern, so an idea can be translate the SP500 on the left,

Interesting isn’t it ? Suppose that we want to forecast (or forsee ?) the SP500 in euro for the next 10 years…

People who enjoy charts would have here a nice tool…

Those two series are extremely correlated, with a correlation of 0.9572,

```> X1=series5[2501:n,4]
> X2=series5[1:(n-2500),5]
> cor(X1,X2)
[1] 0.9572484```

But are the two series cointegrated (see here, here or therefor material on cointegration) ? Well, using standard procedure, we first have to prove that the two series are integrated. First, let us look at the autocorrelograms,

At first sight, we confirm the economic intuition that those indices should be integrated. Standard tests confirm that intuition,

```> acf(X2,lag=1000,col="light green")
> acf(X1,lag=1000,col="light green")
> library(tseries)
Augmented Dickey-Fuller Test
data:  X1
Dickey-Fuller = -1.0768, Lag order = 17, p-value = 0.9264
alternative hypothesis: stationary
Augmented Dickey-Fuller Test
data:  X2
Dickey-Fuller = -1.2905, Lag order = 17, p-value = 0.8788
alternative hypothesis: stationary```

But if we want to go further, we have to find the cointegration relationship between the two series. From an heuristic point of view, a linear regression should be a good proxy,

```> reg=lm(X1~X2)
> plot(residuals(reg))```

`> acf(residuals(reg),lag=1000,col="light green")`

```> adf.test(residuals(reg))
Augmented Dickey-Fuller Test
data:  residuals(reg)
Dickey-Fuller = -5.176, Lag order = 17, p-value = 0.01
alternative hypothesis: stationary
Message d'avis :
In adf.test(residuals(reg)) : p-value smaller than printed p-value
> pp.test(residuals(reg))
Phillips-Perron Unit Root Test
data:  residuals(reg)
Dickey-Fuller Z(alpha) = -46.9775, Truncation lag parameter = 11,
p-value = 0.01
alternative hypothesis: stationary
Message d'avis :
In pp.test(residuals(reg)) : p-value smaller than printed p-value```

When we look at the autocorrelation function, it looks like we do have a stationary series.
This idea is – more or less – the idea of Engle-Granger two step procedure. But actually, we can not directly use Dickey-Fuller’s test to see if residuals are integrated. This was proved in Phillips and Ouliaris (1990), who also proposed a test (see e.g. here),

```> library(tseries); po.test(cbind(X1,X2))
Phillips-Ouliaris Cointegration Test
data:  cbind(X1, X2)
Phillips-Ouliaris demeaned = -53.1766, Truncation lag parameter = 57,
p-value = 0.01
Message d'avis :
In po.test(cbind(X1, X2)) : p-value smaller than printed p-value
Another similar function can be found in R
> library(urca)
> summary(ca.po(cbind(X1,X2)))
########################################
# Phillips and Ouliaris Unit Root Test #
########################################
Test of type Pu
detrending of series none
Call:
lm(formula = z[, 1] ~ z[, -1] - 1)
Value of test-statistic is: 45.2032
Critical values of Pu are:
10pct    5pct    1pct
critical values 20.3933 25.9711 38.3413```

Thus, we has to admit that those series are cointegrated.

Based on that idea, it is possible to model the stationary component, and forecast it for the next ten years, based on the assumption that we know the behavior of one time series. Hence, if we add the confidence interval due to the stationary component uncertainty, we have the following graph,

Of course, again, only uncertainty related to the stationary process is considered here….

# Modélisation des tables de mortalité

En cours d’actuariat, nous commencerons à parler des tables de mortalité prospectives.

Côté théorie, je peux renvoyer vers quelques slides de Pierre Devolder, ici. Je renvoie ici pour un court papier de vulgarisation sur le sujet. Sur l’ajustement de lois de mortalité prospectives, je renvoie au très beau mémoire en ligne ici (ou  pour des choses très proches en anglais). Pour l’approche statique, je renvoie également ici pour une méthode simple d’ajustement de loi de Makeham. Pour les aspects dynamique, il y a des documents techniques ici ou .

Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et  (pour les décès). Il faut retravailler un peu les données, à cause d’une valeur “110+”, mais rien de plus simple,

```> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> DEATH\$Age=as.numeric(as.character(DEATH\$Age))
Message d'avis :
NAs introduits lors de la conversion automatique
> DEATH\$Age[is.na(DEATH\$Age)]=110
> EXPOSURE\$Age=as.numeric(as.character(EXPOSURE\$Age))
Message d'avis :
NAs introduits lors de la conversion automatique
> EXPOSURE\$Age[is.na(EXPOSURE\$Age)]=110```

Coté code, je renvoie à un précédant billet sur le sujet, ici, et  pour le package (non officiel) de démographie développé par Rob Hyndman.

Enfin, à la fin du cours, on parlera un peu de modèles de durée. J’essayerai de taper quelques choses dans les jours à venir, mais côté TD, on utilisera ce que fait Frédéric Planchet (qui connaît fort bien le sujet), et qui se trouve en ligne ici. Les données sont en ligne ici, avec un peu de code R pour retraiter les donnéesici et , et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier ().

Les codes R utilisés en TD sont les suivants (avec en plus les graphiques obtenus),

```> D=DEATH[DEATH\$Year==an,]
> E=EXPOSURE[EXPOSURE\$Year==an,]
> MU = D[,3:5]/E[,3:5]
> plot(0:110,log(MU[,1]),type="l",col="red")
> lines(0:110,log(MU[,2]),col="blue")```

Avec ci-dessus les log des taux de décès instantannés pour l’année 1986.

```> PH=matrix(NA,111,111) # ligne X, colonne H
> PF=PH
> for(x in 0:110){
+ PH[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,2]))
+ PF[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,1]))
+ }
> x=0
> plot(1:111,PH[x+1,],ylim=c(0,1),type="l",col="blue")
> lines(1:111,PF[x+1,],col="red")```

pour la fonction de survie à la naissance (ci-dessus) ou à 50 ans (ci-dessous),

```> somme=function(x){sum(x,na.rm=TRUE)}
> EH=EF=rep(NA,111)
> EH=apply(PH,1,somme)
> EF=apply(PF,1,somme)
> plot(0:110,EH,type="l",col="blue")
> lines(0:110,EF,col="red")```

Ce qui correspond aux espérances de vie résiduelles.

```> for(k in 1:3){
+ Z=MU[,k]
+ Z[is.nan(Z)==TRUE]=NA
+ MU[,k]=Z
+ }
> AGE=unique(DEATH\$Age)
> ANNEE=unique(DEATH\$Year)
> MUF=matrix(MU[,1],length(AGE),length(ANNEE))
> MUH=matrix(MU[,2],length(AGE),length(ANNEE))
> persp(AGE[1:100],ANNEE,log(MUH[1:100,]),
> persp(AGE[1:100],ANNEE,log(MUH[1:100,107:1]),
> image(AGE,ANNEE,log(MUH))```

On a retourné le graph en mettant les années à l’envers, mais on peut aussi regarder la surface pour les femmes, ainsi que l’allure des courbes de niveau

```> persp(AGE[1:100],ANNEE,log(MUF[1:100,107:1]),
> image(AGE,ANNEE,log(MUF))```

Mais nous n’avons ici qu’une lecture statique des tables, et nous ne suivons pas vraiment un individu: en même temps que la personne vieilli, le temps passe…. Considérons une personne née en 1900, et regardons (ex post) l’évolution de sa probabilité de décès (en regardant non pas la probabilité de décéder à 18 d’une personne de 18 ans en 1900, mais en tenant compte du fait que la personne a fêté ses 18 ans en 1918).

```> MUFP=MUHP=rep(NA,111)    # MU PROSPECTIF
> for(i in 1:111){
+     MUHP[i]=MUH[i,colonne+i]
+     MUFP[i]=MUF[i,colonne+i]
+ }
> plot(0:110,log(MUHP),col="purple",type="l")
> D=DEATH[DEATH\$Year==annee.naissance,]
> E=EXPOSURE[EXPOSURE\$Year==annee.naissance,]
> MU = D[,3:5]/E[,3:5]
> lines(0:110,log(MU[,2]),col="blue")```

On peut aussi l’ajustement par un modèle à la Lee-Carter (via la fonction fournie par LifeMetrics)

```> source("https://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
> DEATH=DEATH[DEATH\$Age<90,] # virer age>90
> EXPOSURE=EXPOSURE[EXPOSURE\$Age<90,] # virer age>90
> XV=unique(DEATH\$Age)
> YV=unique(DEATH\$Year)
> ETF=t(matrix(EXPOSURE[,3],length(XV),length(YV)))
> DTF=t(matrix(DEATH[,3],length(XV),length(YV)))
> ETH=t(matrix(EXPOSURE[,4],length(XV),length(YV)))
> DTH=t(matrix(DEATH[,4],length(XV),length(YV)))
> WA=matrix(1,length(YV),length(XV))
> LCF = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)
-1166100 -> -263940.2  ->-204175.2
-188375.9 -> -185822.4  ->-184403.5
-183223.3 -> -183069.8  ->-182350.3
-181701.1 -> -181600.6  ->-181145.5
-180729.6 -> -180671.5  ->-180371.8
-180101.9 -> -180068.6  ->-179871.2
> LCH = fit701(xv=XV,yv=YV,etx=ETH,dtx=DTH,wa=WA)
-2564403 -> -1028776  ->-953629.9
-881334.2 -> -868279.3  ->-834413.6
-804946.8 -> -802817.2  ->-786343
-772775 -> -772007.5  ->-764501.7
-758653.8 -> -758399  ->-755212
-752931 -> -752851  ->-751620.2
-750811.1 -> -750786.7  ->-750351.4
> plot(LCF\$x,LCF\$beta1,type="l",col="red")
> lines(LCH\$x,LCH\$beta1,col="blue")```

```> plot(LCF\$x,LCF\$beta2,type="l",col="red")
> lines(LCH\$x,LCH\$beta2,col="blue")```

```> plot(LCF\$y,LCF\$kappa2,type="l",col="red")
> lines(LCH\$y,LCH\$kappa2,col="blue")```

```> MUF      = DTF/ETF
> MUH      = DTH/ETH
> MUFMODEL = MUF; MUHMODEL = MUH
> for(i in 1:nrow(MUF)){
+ for(j in 1:ncol(MUF)){
+ MUFMODEL[i,j] = exp(LCF\$beta1[j]
+ +LCF\$beta2[j]*LCF\$kappa2[i])
+ MUHMODEL[i,j] = exp(LCH\$beta1[j]
+ +LCH\$beta2[j]*LCH\$kappa2[i])
+ }}```

soit, pour avoir un dessin,

```> persp(LCH\$y,LCH\$x,log(MUH),
> persp(LCH\$y,LCH\$x,log(MUHMODEL),

On peut ainsi comparer les taux bruts (en vert) avec les taux estimés (en bleu)

mais aussi regarder  se trouvent les erreurs le plus grandes,

# Dis moi comment tu t’appelles…

A chaque fois que je prends le métro à Paris, ou que je traverse la gare Montparnasse en descendant du train, je suis surpris par les affiches de publicités qui nous harcèlent. Mais ce qui me surprendra le plus, ce sont les compagnes où on nous présente un Stéphane, un Mathieu, voire un Jean-Pierre, qui doivent nous rendre les affiches plus “humaines” (enfin, je pense). Les gens qui font du marketing ne sont pas forcément stupides (même si personnellement, me montrer la gentille bouille de Stéphane ne va pas me faire aimer quelqu’un qui verbalise des usagers à deux-roues), et je me suis demandé si le choix du prénom (oui, je suis persuadé que le brave Stéphane de l’affiche ne s’appelle pas Stéphane… depuis quand le monde de la pub’ coïnciderait avec le monde réel ? même mes enfants savent que ce n’est pas le cas) était vraiment pertinent…. Bref, j’ai demandé au réseau Quetelet un accès à la base des prénoms…

Pour mon nouvel ami Stéphane, on note effectivement que les Stéphane ont, en moyenne, un peu plus de 36 ans (comme sur la photo), et si on regarde les quantiles, la moitié ont entre 32 ans (quantile à 25%) et 41 ans (quantile à 75%). Bref, on cible bien les trentennaires.

Maintenant, on  peut se demander un peu plus précisément qui est ciblé par cette campagne, on peut regarder d’où viennent les Stéphane,

On retrouve bien un flash pendant une petite dizaine d’années, plutôt dans la partie nord de la France… Bref, la campagne semble très ciblée….

Pour les Mathieu, on veut nous vendre des téléphones portables, je crois…Et effectivement, on est sur un âge moyen beaucoup plus faible, de l’ordre de 25 ans. En revanche, un peu plus étalé que tout à l’heure, avec la moitié ayant entre 15 et 28 ans.

Ici la distribution est plus étalée sur la France, dans le temps,

Bon, le dernier, c’est Jean-Pierre…. alors Jean-Pierre, on ne peut pas trop le voir. On voit qu’il captive des enfants1 en leur racontant son expérience, ce qui correspondrait à une personne un peu âgée, avec de l’expérience. Alors effectivement, si on regarde dans l’historique, les Jean-Pierre ont 55 ans en moyenne, avec 50% entre 50 et 62 ans.

Et là aussi, on retrouve une population proche des Mathieu quant à sa répartition géographique (des pays de cheminots ?).

1 au delà de la fascination exercée par notre ami Jean-Pierre sur ces élèves, pour ma part, ce sont les élèves qui me fascinent… une classe presque exclusivement constituée de gauchers ! peut être pour accentuer le fait qu’il ne faut pas croire que la pub’, c’est comme dans la vraie vie, je ne sais pas… Remarquez, je n’ai rien contre les gauchers sur les photos ! mais c’est un peu comme les auvergnats: un ça va, mais c’est quand il y en a trop que ça devient gênant….

# Vendredi 13 et probabilités….

Je ne sais pas si c’est la proximité avec Halloween, mais je trouvais ce vendredi 13 qui s’annonce avait quelque chose de spécial.. Et tant pis si après mes études sur les jours de naissance les jours de pleine lune, je passe définitivement pour un mystique. Disons plutôt que j’aime bien interroger les croyances de grands mères…. de la même manière que j’aime questionner les pratiques des actuaires (souvent pleines de bon sens, mais parfois un peu bancales).
Le problème avec les vendredi 13 est que l’on ne sait jamais si c’est un jour qui porte chance, ou au contraire un jour qui doit porter malheur… Dans le doute, je n’ai pas pu m’empêcher de mettre quelques trèfles à quatre feuilles, des brin de muguet et un pompon de marin.. Comme disait un très grand scientifique, “je n’y crois pas, mais il parait que ça marche même pour ceux qui n’y croient pas“. Je ne vais donc pas me priver !

Bon notons tout d’abord qu’avoir un vendredi 13 n’a rien d’exceptionnel… il y en a déjà au moins un par an. Mais ce qui est troublant c’est que le 13 d’un mois a plus de chances de tomber un vendredi qu’un autre jour. Pour avoir un cycle de semaines entières, il faut compter 400 ans, soit 20871 semaines (ce que l’on appelle un cycle grégorien). Mais en 400 ans, on a 4800 mois, soit 4800 “13”… qui n’est pas divisible par 7, ce qui tend à laisser penser que certains jours de la semaine tomberont (un peu) plus souvent que les autres… Et Brown a montré que c’était le cas (dès 1933, dans le American Mathemical Monthly1). Pour les amateurs, Jean-Luc Nothias avait fait un article en mars dernier sur le sujet (ici).
De plus, certains auteurs ont tenté des études statistiques sur le sujet. Par exemple, en 1993, le British Medical Journal avait publié un article sur le sujet (“Is Friday the 13th Bad for Your Health?” , ici) “Friday 13th is unlucky for some. The risk of hospital admission as a result of a transport accident may be increased by as much as 52 percent. Staying at home is recommended.” Bon, on va essayer de jeter un oeil….

• Visites chez le médecin (généraliste) les vendredi 13

Les données a notre disposition (gentillement données par la Caisse d’Assurance Maladie) correspondent au nombre de consultations d’un médecin généraliste, par jour… Bon, ces données présentent un soucis technique: les données arrivent avec un (léger) délai, et les données du vendredi arrivent regroupées avec celles du samedi et du dimanche le mardi suivant. On ne peut donc pas comparer un vendredi au jeudi précédant par exemple… Pour éviter ce genre de soucis, la stratégie la plus simple consiste à comparer le nombre de consultations chez le médecin….

avec en rouge les vendredi 13, et en jaune la moyenne des vendredi 6 et 20…. Vu de loin, on consulte un peu plus les vendredi 13 que les vendredi avant et après. Si on fait un test de Student (quelle que soit l’hypothèse alternative) on rejette l’hypothèse comme quoi on conculte plus son médecin un vendredi 13…

```> t.test(V13,V)
Welch Two Sample t-test
data:  V13 and V
t = 0.6126, df = 9.052, p-value = 0.5552
alternative hypothesis: true difference in means is not equal to 0```

Mais peut être est-ce par qu’une consultation chez le médecin généraliste n’est pas associé à de la “malchance” ?

• Quid des catastrophes naturelles ?

Cherchons autre chose que l’on pourrait associer à cette prétendue “malchance“. Si on regarde les communes classées catastrophes naturelles en France, entre août 1982 et mars 2008, il y a eu 131 285 classements de communes comme touchées par une catastrophe naturelle.(et paru au Journal Officiel, ces données sont publiques). Seuls 393 classements ont été  fait un vendredi 13, ce qui pourrait sembler très faible… En fait, cette série est très ératique.

Le graphique montre l’évolution du nombre de classements les mois où il y a eu un vendredi 13. Statistiquement, on pourrait noter qu’il y a beaucoup plus de communes touchées par une catastrophe naturelle un mercredi 11 par exemple… Bref, la base semble inutilisable pour cette étude… simplement car un cinquième des classements concernent les tempêtes de décembre 1999 (ou plutôt les inondations qui ont suivies, car les tempêtes sont exclues du régime cat nat, en métropole).

• Quid des gros sinistres de perte d’exploitation et des incendies entreprises…

J’aurais préféré avoir d’autres données d’assureurs, mais comme toujours, lorsque j’ai demandé les données, on m’a expliqué que c’étaient des données sensibles, et que je ne pourrais rien avoir ! Fort heureusement, il me reste quelques jeux de données…. en particulier une base de gros sinistres de sinistres en perte d’exploitation en France.
Cette fois, ma base va de début 1985 à fin 2000, avec en tout 2387 sinistres (cette fois l’interpréation devrait être plus simple car les données sont identifiées comme des sinistres). Au final 17 sinistres majeurs2 se sont produits un vendredi 13. Cette fois, le chiffre est significativement grand, car par exemple, seulement 4 se sont produit un jeudi 12, et 9 un samedi 14…

Là encore, si on regarde ce qui s’est passé les mois où il y avait un vendredi 13, on note qu’il y a eu beaucoup de sinistres ces jours là… Envie de faire un test ? On a en tout 2387 sinistres à faire tomber sur 5824 jours, parmi lesquels, on compte 27 vendredi 13. Bref, en supposant les sinistres indépendants, chaque sinistre à une probabilité de 0,463% de tomber un vendredi 13. Le nombre de sinistres tombant un vendredi 13 est donc supposé suivre une loi binomiale

On pourrait espérer avoir 11 sinistres qui tombent un vendredi 13 (et on en observe 17). Si on calcule le quantile à 95%, on tombe sur 17. En fait, on avait environ 1 chance sur 20 d’avoir autant de sinistres ces jours là…. ce qui correspond exactement au seuil que l’on retient pour juger de la significativité d’un test (comme je l’évoquais un peu ici). Bref, on ne peut pas vraiment conclure….
En fait, comme il s’agissait de sinistres d’entreprises, je pensais qu’il y aurais un biais, car (par définition) les vendredi 13 tombent un vendredi, et que les jours de semaines devaient avoir plus de sinistres que les jours de week end. Mais ce n’est pas le cas

Bref, on est dans la zone d’incertitude des tests statistiques….
Mais on peut aussi utiliser une base d’incendies intreprises que j’ai sous le coude depuis quelques années, avec des sinistres survenus entre 1982 et 1996. Cette base contient 9613 gros sinistres: En utilisant la même méthode que pour les sinistres de perte d’exploitation, on note que beaucoup d’incendies sont survenus un vendredi 13 dans cette base: 52. C’est un peu plus que ce que nous aurions pu espérer, à savoir 45,62. Là encore, en notant que le nombre d’incendies devrait suivre une loi binomiale, notons qu’il y avait 15,4% de chances d’atteindre 52 sinistres, ce qui n’est – cette fois-ci – pas particulièrement rare.

Par contre, sur cette base, l’effet week end semble plus prononcé, ce qui semble aller dans le sens que ces 52 sinistres, c’est loin d’être exceptionnellement elévé pour un jour de semaine.

On pourrait aussi parler de karma, ou de chose dans le genre….

1 B. H. Brown, “Solution to Problem E36”, American Mathematical Monthly, vol. 40, issue 10 (1933), p. 607 Pour faire simple, un cycle grégorien correspond à 400 années, soit 136097 jours, ou encore 20 871 semaines. Si on regarde la distribution des jours on voit que samedi et jeudi tombent 684 fois, lund et mardi 685, mercredi et dimanche 687, et vendredi 688. Bon, d’accord, je pinaille…

2 pour les plus curieux, ou les plus sceptiques, je mets ici les codes pour ceux qui veulent jouer avec les données… c’est un exercice amusant pour jouer avec la manipulation de dates sous R…

```> B1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/data-cat-nat.txt",header=TRUE)
> B2\$DATEPE=as.character(B2\$DATEPE)
> x2=as.Date(B2\$DATEPE,"%d/%m/%Y")
> B3\$date=as.character(B2\$date)
> x3=as.Date(B3\$date,"%d/%m/%Y")
> p2=sum((x3>=min(x2))&(x3<=max(x2)))/as.numeric(max(x2)-min(x2))
> n2=length(x2);  n2*p2
> qbinom(.05,n2,p2); qbinom(.95,n2,p2)
> pbinom(sum(x2%in%x3)  ,n2,p2)
> NOM=weekdays(x3[1]+1:7-5)
> S=rep(NA,7)
> for(i in 1:7){S[i]=sum(weekdays(x2)==NOM[i])}
> barplot(S,names.arg=NOM,col="light blue")```

voilà pour les sinistres en perte d’exploitation… le reste n’est qu’une variation sur le même thème…

# Séminaire Probabilité et Statistique, UBO, Brest

Talk at the statistical seminar at the Université de Bretagne Occidentale, in Brest, Wednesday May 6th Tuesday May 5th, 14h (in  10 days), on “multivariate extremes. Slides can be found here.

The talk will give a detailed introduction on multivariate extremes and related concepts. Then the case of Archimedean copula will be fully described (following the paper with Johan Segers).

[04/05/2009]: some applications in risk management will be shown at the end of talk, as well as some news things on spatial correlation.

and in order to illustrate tail convergence of Archimedean copulas, I have uploaded two animations, with tail independence below,

with tail dependence (or asymptotic dependence),

# Short course on catastrophes: the Mexican earthquake

One case study of catastrophes covers will be the Mexican cat bond. The dataset of past historical events in Mexico can be downloaded here.

Since everything we need is in this dataset, everything can be done using R.

```mex=read.table("https://blogperso.univ-rennes1.fr/arthur.charpentier/public/conferences/EQ-mex-6-5.txt")
mex\$D=mex\$V1*10000+mex\$V2*100+mex\$V3
mex\$Date=as.Date(as.character(mex\$D),"%Y%m%d")```

# Convolution(s) et lois composées, aspects numériques

Nous avons vu en cours que l’outil de base pour travailler sur les sommes de variables aléatoires était la convolution,

 $(f * g )(t) ,$ $stackrel{mathrm{def}}{=} int_{-infty}^{infty} f(tau)cdot g(t - tau), dtau$ $= int_{-infty}^{infty} f(t-tau)cdot g(tau), dtau.$

permettant d’obtenir la densité ou la fonction de répartition de la somme de deux variables indépendantes.

Cet outils permet en théorie d’obtenir la forme d’une loi composée, i.e. d’une somme dont le nombre de termes est aléatoire.

Pour les calculs numériques de la convolée, Vincent Goulet a programmé pas  mal de choses dans la librairie actuar. Numériquement, il faut discrétiser la loi des coûts de sinistres. Si on suppose par exemple que les coûts individuels sont discrets, on peut obtenir la distribution simplement

```> fx <- c(0, 0.15, 0.2, 0.25, 0.125, 0.075,
+         0.05, 0.05, 0.05, 0.025, 0.025)
> pn <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.15, 0.06, 0.03, 0.01)
> Fs <- aggregateDist("convolution", model.freq = pn,
+                     model.sev = fx, x.scale = 25)
> summary(Fs)
Aggregate Claim Amount Empirical CDF:
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.0   150.0   275.0   314.5   400.0  2000.0
> c(Fs(0), diff(Fs(25 * 0:21))) # probability mass function
[1] 0.05000000 0.01500000 0.02337500 0.03467500 0.03257656 0.03578639
[7] 0.03980787 0.04356232 0.04751800 0.04903380 0.05189806 0.05137886
[13] 0.05118691 0.05030486 0.04818189 0.04575882 0.04280890 0.03937836
[19] 0.03574568 0.03196808 0.02832446 0.02478833
> plot(Fs)```

On peut aussi utiliser des méthodes d’approximation récursives,

> Fs <- aggregateDist(“recursive”, model.freq = “poisson”,
+                     model.sev = fx, lambda = 3, x.scale = 25)
> plot(Fs)

Ce qui donne le graphique suivant

Des méthodes d’approximation de type normal power sont aussi programmées, ainsi que les techniques de simulations. Par exemple ici pour une survenance Poissonienne avec des coûts Gamma, on écrit

```> Fs <- aggregateDist("recursive", model.freq = "poisson",
+ model.sev = fx, lambda = 3, x.scale = 25)
> plot(Fs)```

Enfin, pour lier avec le DM que je donnerais mercredi, Vincent Goulet a aussi programmé une méthode d’estimation de quantile,

```> model.freq <- expression(data = rpois(3))
> model.sev <- expression(data = rgamma(100, 2))
> Fs <- aggregateDist("simulation", nb.simul = 1000,
+                     model.freq, model.sev)```

# Estimation de quantile par noyau beta

Le papier sur l’estimation de quantile par noyau beta, coécrit avec Abder Oulidi, est accepté pour publication dans Statistics and Computing, http://link.springer.com/…

In this paper we propose several nonparametric estimators of quantiles based on Beta kernel and applied to transformed data by the generalized Champernowne distribution initially fitted to the data. A Monte-Carlo based study will show that those estimators improve the efficiency of a traditional ones, not only for light tailed distributions, but also heavy tails, when the probability level is close to 1. We also compare these estimators with the Extreme Value Theory Quantile applying to Danish data on large fire insurance losses.