Promesse et rationalité économique (à court terme)

Non, je vais pas faire une analyse des promesses en cette année électorale (bien que le sujet pourrait être amusant, on peut penser aux propos du candidat Newt Gingrich qui veut que la lune devienne un état américain, ici ou ). Non, en fait je voulais faire un billet sur le sujet les assureurs sont-ils vraiment des voleurs ? En effet, lorsqu’en classe on présente le fonctionnement d’un contrat d’assurance – en particulier lorsque l’on parle de tarification et de l’inversion du cycle de production, ou du provisionnement – on explique que l’assuré achète une promesse: celle de se faire indemniser les sinistres survenus pendant la période dite de couverture. Car formellement, c’est bien de cela qu’il s’agit. Sauf qu’il y a toujours quelqu’un pour faire noter que ce n’est peut-être pas aussi simple et que l’assureur, souvent, rechigne à payer.

En fait, il est parfaitement rationnel – d’un point de vue économique (on s’entend) – que l’assureur négocie. En tous les cas si l’assureur a une vision à très court terme… Et c’est ce que justifient Briegleb & Lemaire (1982), sur l’analyse du marchandage ou Lemaire (1982) (dont sera très largement inspiré ce billet).
Considérons un cas très simple: un assuré a subit une perte http://freakonometrics.blog.free.fr/public/perso5/tribunal33.gif (connue, de l’assuré comme de l’assureur). Mais l’assureur refuse de payer. Et l’assuré envisage d’aller au tribunal. Afin de voir s’il s’agit d’une décision rationnelle, on suppose que l’assuré (http://freakonometrics.blog.free.fr/public/perso5/tribunal18.gif) avait une richesse  avant que le sinistre ne survienne, qu’aller en justice lui coûte , et qu’il estime que sa probabilité de gagner le procès est http://freakonometrics.blog.free.fr/public/perso5/tribunal05.gif. De son coté, la compagnie (http://freakonometrics.blog.free.fr/public/perso5/tribunal19.gif) d’assurance a une richesse , qu’aller en justice lui coute http://freakonometrics.blog.free.fr/public/perso5/tribunal04.gif, et que la probabilité que la compagnie perde le procès est – selon la compagnie – http://freakonometrics.blog.free.fr/public/perso5/tribunal06.gif (les http://freakonometrics.blog.free.fr/public/perso5/tribunal.gif sont alors les probabilités que l’assuré gagne son procès, mais vu des deux cotés). Au tribunal, on suppose que le juge rend un verdict du genre tout ou rien (on n’autorise pas ici de compensation plus importantes que le cout réel du sinistre, et on ne parle pas non plus de remboursement des frais d’avocat en cas de victoire du procès).
L’espérance d’utilité de l’assuré est

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

alors que pour la compagnie d’assurance, son espérance d’utilité est

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

La situation est Pareto optimale si

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

Si on suppose que les agents sont risquophobe (i.e. http://freakonometrics.blog.free.fr/public/perso5/tribunal10.gifdécroissante), alors

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

alors que

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

Aussi, une condition nécessaire pour que la situation soit Pareto-optimale est que

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

i.e. http://freakonometrics.blog.free.fr/public/perso5/tribuanl16.gif. Aller au tribunal n’est pas Pareto optimal si http://freakonometrics.blog.free.fr/public/perso5/tribunal18.gif croit plus en ses raisons de gagner que http://freakonometrics.blog.free.fr/public/perso5/tribunal19.gif. Si http://freakonometrics.blog.free.fr/public/perso5/tribunal17.gif, il peut alors être optimal de ne pas aller au tribunal…. On ne dit pas ici qu’il faut que la compagnie indemnise (intégralement) l’assuré, mais que l’assureur et l’assuré peuvent avoir intérêt a trouver une solution amiable…
Si http://freakonometrics.blog.free.fr/public/perso5/tribuanl16.gif, il peut être intéressant d’aller au tribunal… mais pas forcément… Du point de vue de l’assureur, on peut supposer que sa richesse lui permet de ne pas être risquophobe, ou au moins d’être risque neutre (face a cette assuré au moins). Elle est donc indifférente entre aller au tribunal (et avoir un rendement incertain) et verser un règlement amiable si le versement correspond à son espérance de perte si elle va au procès. Autrement dit http://freakonometrics.blog.free.fr/public/perso5/tribunal20.gif. L’assuré acceptera cette somme si l’utilité qu’il en tire est plus grande que l’espérance d’utilité qu’il aura en allant au tribunal,

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

La différence entre le terme de gauche et le terme de droite est une fonction http://freakonometrics.blog.free.fr/public/perso5/tribunal22.gif. Si l’agent a une utilité exponentielle

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

alors

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

doit être une fonction positive. Un calcul rapide de dérivées montre que cette fonction est convexe, avec un minimum atteint en

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

Moralité ?

  • pour les petits et pour les gros sinistres, il est normal que l’assureur et l’assuré cherchent une solution amiable (transaction)
  • en revanche, pour les sinistres de cout intermédiaire, il peut être légitime d’aller au tribunal (procès).

Maintenant, si on a les grandes directions… peut-être peut-on essayer de mieux comprendre quelle serait le montant de la transaction dans le cas ou une solution amiable est envisagée. Nash (1950) a envisagé un jeu non-coopératif afin de décrire le marchandage. La transaction optimale est la quantité http://freakonometrics.blog.free.fr/public/perso5/tribunal28.gif qui maximise le produit des deux gains d’utilité (cf Nash (1950) ou Roth (1979))

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

La condition du premier ordre donne l’équation suivante

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

Si on fait un peu de calcul numérique, on note que la solution http://freakonometrics.blog.free.fr/public/perso5/tribunal32.gif est décroissante en http://freakonometrics.blog.free.fr/public/perso5/tribunal31.gif: plus l’assuré est risquophobe, plus faible sera la montant de la transaction…. Ce qui n’est pas surprenant, on retrouve ici un résultat noté dans Peeters (1981), Roth (1982) ou ecnore Osbourne (1985). Maintenant, sur le modèle théorique on pourrait bien entendu argumenter que le tribunal fonctionne selon une règle trop simple du tout ou rien, et qu’un procès peut couter beaucoup cher que le simple remboursement du cout du sinistre. Ou qu’il pourrait créer une jurisprudence. On pourrait aussi regarder économétriquement le lien entre le cout des sinistres et les montants des ententes amiables (c’est ce que fait Ayusoa, Bermúdezb & Santolinoc (2012)) mais c’est une autre histoire…. ou ça pourrait être le sujet pour une autre histoire… à suivre donc.

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J’avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu’on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n’est pas forcément l’unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c’est d’ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), “when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.”

N’étant toujours pas convaincu, j’évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

http://freakonometrics.hypotheses.org/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")

http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l’impact de http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date http://freakonometrics.hypotheses.org/files/2015/12/lex14.gif, pour un horizon http://freakonometrics.hypotheses.org/files/2015/12/lisse12.gif est alors http://freakonometrics.hypotheses.org/files/2015/12/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1

http://freakonometrics.hypotheses.org/files/2015/12/lisse14.gif

> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c’est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.

Call:
HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])

Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE

Coefficients:
[,1]
a 805.0389

> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.

http://freakonometrics.hypotheses.org/files/2015/12/lex50.gif

et

http://freakonometrics.hypotheses.org/files/2015/12/HW02.gif

Dans ce cas, la prédiction est http://freakonometrics.hypotheses.org/files/2015/12/lisse13.gif. Le code pour faire un lissage double est là encore assez simple,

> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }

Sur la figure suivante, on visualise l’évolution du lissage en fonction de http://freakonometrics.hypotheses.org/files/2015/12/lisse15.gif (http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif étant ici fixé),
http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-2.gif
(le lissage simple – avec le même poids http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Visualization in regression analysis

Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5),

> library(gdata)
> XLS1=read.xls("http://api.worldbank.org/datafiles
/NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")]
> names(data1)[3]="GDP"
> XLS2=read.xls("http://api.worldbank.org/datafiles
/SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data2=XLS2[-(1:28),c("Country.Code","X2010")]
> names(data2)[2]="MORTALITY"
> data=merge(data1,data2)
> head(data)
Country.Code         Country.Name       GDP MORTALITY
1          ABW                Aruba        NA        NA
2          AFG          Afghanistan  1207.278     149.2
3          AGO               Angola  6119.930     160.5
4          ALB              Albania  8817.009      18.4
5          AND              Andorra        NA       3.8
6          ARE United Arab Emirates 47215.315       7.1

If we estimate a simple linear regression – http://freakonometrics.blog.free.fr/public/perso5/logormal01.gif  – we get

> regBB=lm(MORTALITY~GDP,data=data)
> summary(regBB)

Call:
lm(formula = MORTALITY ~ GDP, data = data)

Residuals:
Min     1Q Median     3Q    Max
-45.24 -29.58 -12.12  16.19 115.83

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.1008781  4.1577411  16.139  < 2e-16 ***
GDP         -0.0017887  0.0002161  -8.278 3.83e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 39.99 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.2909,	Adjusted R-squared: 0.2867
F-statistic: 68.53 on 1 and 167 DF,  p-value: 3.834e-14

We can look at the scatter plot, including the linear regression line, and some confidence bounds,

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5)",cex=.5)
> text(data$GDP,data$MORTALITY,data$Country.Name,pos=3)
> x=seq(-10000,100000,length=101)
> y=predict(regBB,newdata=data.frame(GDP=x),
+ interval="prediction",level = 0.9)
> lines(x,y[,1],col="red")
> lines(x,y[,2],col="red",lty=2)
> lines(x,y[,3],col="red",lty=2)

We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,

> boxcox(regBB)

it looks like taking the logarithm of the mortality rate should be better, i.e. http://freakonometrics.blog.free.fr/public/perso5/lognormal02.gif or http://freakonometrics.blog.free.fr/public/perso5/lognormal05.gif:

> regLB=lm(log(MORTALITY)~GDP,data=data)
> summary(regLB)

Call:
lm(formula = log(MORTALITY) ~ GDP, data = data)

Residuals:
Min      1Q  Median      3Q     Max
-1.3035 -0.5837 -0.1138  0.5597  3.0583

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.989e+00  7.970e-02   50.05   <2e-16 ***
GDP         -6.487e-05  4.142e-06  -15.66   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7666 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.5949,	Adjusted R-squared: 0.5925
F-statistic: 245.3 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5,log="y")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=seq(300,100000,length=101)
> y=exp(predict(regLB,newdata=data.frame(GDP=x)))*
+ exp(summary(regLB)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scale or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5)

on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.

But why shouldn’t we take also the logarithm of the GDP ? We can fit a model http://freakonometrics.blog.free.fr/public/perso5/lognormal03.gif or equivalently http://freakonometrics.blog.free.fr/public/perso5/lognormal04.gif.

> regLL=lm(log(MORTALITY)~log(GDP),data=data)
> summary(regLL)

Call:
lm(formula = log(MORTALITY) ~ log(GDP), data = data)

Residuals:
Min       1Q   Median       3Q      Max
-1.13200 -0.38326 -0.07127  0.26610  3.02212

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.50192    0.31556   33.28   <2e-16 ***
log(GDP)    -0.83496    0.03548  -23.54   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5797 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.7684,	Adjusted R-squared: 0.767
F-statistic:   554 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5,log="xy")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=exp(seq(1,12,by=.1))
> y=exp(predict(regLL,newdata=data.frame(GDP=x)))*
+ exp(summary(regLL)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scales or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5)

on the standard scale. If we compare the last two predictions, we have

with in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).

Examen intra, modèles de régression

L’examen intra comportera deux exercices proches de ceux faits en démonstration, et une analyse de sorties de régression. Cette dernière partie portera sur la base suivante (ou peut être une extraction de cette base)

> examen=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv",
+ header=TRUE,sep=";")
> tail(examen)
EQUIPE halfdiff winner
9917                 MARYLAND       -1      0
9918             RHODE ISLAND      -15      1
9919            ROBERT MORRIS        4      1
9920 NORTH CAROLINA-ASHEVILLE       11      1
9921                    BROWN       -2      1
9922          MOUNT ST. MARYS        1      0
totalpointsdiff secondhalfpoints
9917             -10               26
9918               1               44
9919               6               43
9920              19               39
9921               6               48
9922              -9               26
nodiff_pts_1000_left hometeam halfwinner
9917                    9        1          0
9918                   18        1          0
9919                   17        1          1
9920                   20        1          1
9921                   16        1          0
9922                   13        1          1
halflooser secondhalfpointsother
9917          1                    35
9918          1                    28
9919          0                    41
9920          0                    31
9921          1                    40
9922          0                    36

qui contient des matchs de baskets universitaires.

  • EQUIPE le nom de l’équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
  • WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps

Tests on tail index for extremes

Since several students got the intuition that natural catastrophes might be non-insurable (underlying distributions with infinite mean), I will post some comments on testing procedures for extreme value models.

A natural idea is to use a likelihood ratio test (for composite hypotheses). Let http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif denote the parameter (of our parametric model, e.g. the tail index), and we would like to know whether http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif is smaller or larger than http://freakonometrics.blog.free.fr/public/perso5/lrtest22.gif (where in the context of finite versus infinite mean http://freakonometrics.blog.free.fr/public/perso5/lrtest23.gif). I.e. either http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif belongs to the set http://freakonometrics.blog.free.fr/public/perso5/lrtest-10.gif or to its complementary http://freakonometrics.blog.free.fr/public/perso5/lrtest-11.gif. Consider the maximum likelihood estimator http://freakonometrics.blog.free.fr/public/perso5/lrtest24.gif, i.e.

http://freakonometrics.blog.free.fr/public/perso5/lrtest-9.gif

Let http://freakonometrics.blog.free.fr/public/perso5/lrtest25.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-3.gif denote the constrained maximum likelihood estimators on http://freakonometrics.blog.free.fr/public/perso5/lrtest26.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest27.gif respectively,

http://freakonometrics.blog.free.fr/public/perso5/lrtest-12.gif

http://freakonometrics.blog.free.fr/public/perso5/lrtest-2.gif

Either http://freakonometrics.blog.free.fr/public/perso5/lrtest-13.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-6.gif (on the left), or http://freakonometrics.blog.free.fr/public/perso5/lrtest-14.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-7.gif (on the right)

So likelihood ratios

http://freakonometrics.blog.free.fr/public/perso5/lrtest-15.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-16.gif

 are either equal to

http://freakonometrics.blog.free.fr/public/perso5/lrtest-19.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-18.gif

or

http://freakonometrics.blog.free.fr/public/perso5/lrtest-20.gif        http://freakonometrics.blog.free.fr/public/perso5/lrtest-17.gif

If we use the code mentioned in the post on profile likelihood, it is easy to derive that ratio. The following graph is the evolution of that ratio, based on a GPD assumption, for different thresholds,

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM
> U=seq(2,10,by=.2)
> LR=P=ES=SES=rep(NA,length(U))
> for(j in 1:length(U)){
+ u=U[j]
+ Y=X[X>u]-u
+ loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
+ XIV=(1:300)/100;L=rep(NA,300)
+ for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
+ plot(XIV,L,type="l")
+ PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
+ (L0=(OPT=optimize(f=PL,interval=c(0,10)))$objective)
+ profilelikelihood=function(beta){
+ -loglikelihood(1,beta) }
+ (L1=optim(par=1,fn=profilelikelihood)$value)
+ LR[j]=L1-L0
+ P[j]=1-pchisq(L1-L0,df=1)
+ G=gpd(X,u)
+ ES[j]=G$par.ests[1]
+ SES[j]=G$par.ses[1]
+ }
>
> plot(U,LR,type="b",ylim=range(c(0,LR)))
> abline(h=qchisq(.95,1),lty=2)

with on top the values of the ratio (the dotted line is the quantile of a chi-square distribution with one degree of freedom) and below the associated p-value

> plot(U,P,type="b",ylim=range(c(0,P)))
> abline(h=.05,lty=2)

In order to compare, it is also possible to look at confidence interval for the tail index of the GPD fit,

> plot(U,ES,type="b",ylim=c(0,1))
> lines(U,ES+1.96*SES,type="h",col="red")
> abline(h=1,lty=2)

To go further, see Falk (1995), Dietrich, de Haan & Hüsler (2002), Hüsler & Li (2006) with the following table, or Neves & Fraga Alves (2008). See also here or there (for the latex based version) for an old paper I wrote on that topic.

Lissage exponentiel et séries temporelles

L’idée du lissage exponentiel est de supposer que http://freakonometrics.hypotheses.org/files/2015/12/lex01.gif, où on suppose que le bruit http://freakonometrics.hypotheses.org/files/2015/12/lex02.gif est ici un bruit indépendant (imprévisible). On va alors supposer que

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

où http://freakonometrics.hypotheses.org/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif.

On parle de lissage pour moyenne mobile à poids exponentiel (ewma, exponentially weighted moving average, évoqué ici ou  sur ce blog). En effet

http://freakonometrics.hypotheses.org/files/2015/12/lex08.gif

i.e., en itérant

http://freakonometrics.hypotheses.org/files/2015/12/lex09.gif

(même s’il faut faire attention à cette décomposition qui, en pratique, s’arrête avec la première observation).

La méthode de Holt-Winters propose de généraliser ce lissage, en introduisant deux (voire trois) composantes.

  • http://freakonometrics.hypotheses.org/files/2015/12/lex05.gif sera interprété comme un niveau
  • http://freakonometrics.hypotheses.org/files/2015/12/lex11.gif correspondant à la variation du processus
  • http://freakonometrics.hypotheses.org/files/2015/12/lex12.gif qui correspond à un cycle (éventuel)

On suppose désormais, dans cette méthode (dite de Holt-Winters), que

http://freakonometrics.hypotheses.org/files/2015/12/lex50.gif

et

http://freakonometrics.hypotheses.org/files/2015/12/HW02.gif

D’un point de vue plus pratique, la programmation se fait comme suit,

> HW=HoltWinters(X,alpha=.2,beta=0)
> plot(HW)

Cette trois composantes seront utilisées pour faire une prévision de la série: la prédiction sera de la forme http://freakonometrics.hypotheses.org/files/2015/12/lex16.gif, à la date http://freakonometrics.hypotheses.org/files/2015/12/lex14.gif, pour un horizon http://freakonometrics.hypotheses.org/files/2015/12/lex17.gif (avec les notations du précédant billet).

> P=predict(HW,24,prediction.interval=TRUE)
> plot(HW,xlim=range(c(time(X),time(P))))
> polygon(c(time(P),rev(time(P))),c(P[,2],rev(P[,3])),
+ col="yellow",border=NA)
> lines(P[,1],col="red",lwd=3)

What is working as an actuary really like?

Just a short post to mention the page entitled “What is working as an actuary really like ?” on http://beanactuary.com/. I wanted to share that page with some students (even if it is quite incomplete) because I am always surprised when I hear that what we teach during lectures is too mathematical (or even to much based on computing skills). Practice can be much more complex than what we do in our courses (e.g. variable selection in ratemaking, integrating expert opinion in claims reserving, or hedging life insurance contracts with several complex financial options). If practitioners want to share their experience, comments are welcome… I just hope that some students will stop thinking that Barney Stinson is doing the job of an actuary… and that what we teach might – somehow – be useful someday…

the Dirichlet distribution

In the course, since we are still introducing some concepts of dependent distributions, we will talk about the Dirichlet distribution, which is a distribution over the simplex of http://freakonometrics.hypotheses.org/files/2017/07/diri11.gif. Let http://freakonometrics.hypotheses.org/files/2017/07/diri01.gif denote the Gamma distribution with density (on http://freakonometrics.hypotheses.org/files/2017/07/diri03.gif)

http://freakonometrics.hypotheses.org/files/2017/07/diri02.gif

Let http://freakonometrics.hypotheses.org/files/2017/07/diri04.gif denote independent http://freakonometrics.hypotheses.org/files/2017/07/diri05.gif random variables, with http://freakonometrics.hypotheses.org/files/2017/07/diri06.gif. Then http://freakonometrics.hypotheses.org/files/2017/07/diri07.gif where

http://freakonometrics.hypotheses.org/files/2017/07/diri08.gif

has a Dirichlet distribution with parameter

http://freakonometrics.hypotheses.org/files/2017/07/diri09.gif

Note that http://freakonometrics.hypotheses.org/files/2017/07/diri10.gif has a distribution in the simplex of http://freakonometrics.hypotheses.org/files/2017/07/diri11.gif,

http://freakonometrics.hypotheses.org/files/2017/07/diri40.gif

and has density

http://freakonometrics.hypotheses.org/files/2017/07/diri12.gif

We will write http://freakonometrics.hypotheses.org/files/2017/07/diri13.gif.

The density for different values of http://freakonometrics.hypotheses.org/files/2017/07/diri20.gif can be visualized below, e.g. http://freakonometrics.hypotheses.org/files/2017/07/diri21.gif, with some kind of symmetry,
http://freakonometrics.hypotheses.org/files/2017/07/dirichlet222.gif
or http://freakonometrics.hypotheses.org/files/2017/07/diri22.gif and http://freakonometrics.hypotheses.org/files/2017/07/diri23.gif, below
http://freakonometrics.hypotheses.org/files/2017/07/dirichlet522.gif
and finally, below, http://freakonometrics.hypotheses.org/files/2017/07/diri24.gif


Note that marginal distributions are also Dirichlet, in the sense that if

http://freakonometrics.hypotheses.org/files/2017/07/diri13.gif

then

http://freakonometrics.hypotheses.org/files/2017/07/diri14.gif

if http://freakonometrics.hypotheses.org/files/2017/07/diri15.gif, and if http://freakonometrics.hypotheses.org/files/2017/07/diri16.gif, then http://freakonometrics.hypotheses.org/files/2017/07/diri17.gif‘s have Beta distributions,

http://freakonometrics.hypotheses.org/files/2017/07/diri18.gif

See Devroye (1986) section XI.4, or Frigyik, Kapila & Gupta (2010) .This distribution might also be called multivariate Beta distribution. In R, this function can be used as follows

> library(MCMCpack)
> alpha=c(2,2,5)
> x=seq(0,1,by=.05)
> vx=rep(x,length(x))
> vy=rep(x,each=length(x))
> vz=1-x-vy
> V=cbind(vx,vy,vz)
> D=ddirichlet(V, alpha)
> persp(x,x,matrix(D,length(x),length(x))

(to plot the density, as figures above). Note that we will come back on that distribution later on so-called Liouville copulas (see also Gupta & Richards (1986)).

Visualisation des résidus pour des données spatiales

Mardi, nous allons travailler un peu sur la modélisation du nombre d’homicides aux États-Unis, à partir de la base

> US=read.table("http://freakonometrics.free.fr/US.txt", 
+ header=TRUE,sep=";")

(je renvoie sur le précédant billet pour un descriptif précis). Idéalement, ça serait parfait si on pouvait visualiser sur une carte les variables. Pour cela, il faut rajouter une colonne à notre base, avec le nom complet des états,

> abreviation=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/etatus.csv", 
> header=TRUE,sep=",") 
> US$USPS=rownames(US) 
> US=merge(US,abreviation) 
> US$nom=tolower(US$NOM)

Cette fois, on va pouvoir faire de la cartographie, les noms de nos états étant (presque) les mêmes que ceux des cartes de R,

> library(maps) 
> VL0=strsplit(map("state")$names,":") 
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,US$nom)

Cette fois-ci, on a toutes les informations pour faire une carte, avec une couleur fonction de la variable d’intérêt (espérance de vie à la naissance, taux d’homicides, taux illettrisme, etc).

> library(RColorBrewer) 
> carte=function(V=US$Murder,titre= 
+ "Taux d'homicides aux Etats-Unis"){ 
+ variable=as.numeric(as.character(cut(V, 
+ quantile(V,seq(0,1,by=1/6)),labels=1:6))) 
+ niveau=variable[ETAT] 
+ couleur=rev(brewer.pal(6,"RdBu")) 
+ noml=levels(cut(V,quantile(V,seq(0,1,by=1/6)))) 
+ map("state", fill = TRUE, col=couleur[niveau]); 
+ legend(-78,34,legend=noml,fill=couleur, + cex=1,bty="n"); 
+ title(titre)}

Commençons par l’analyse du nombre de jours durant lesquels la température passe en dessous de 0°C, par an, en moyenne, dans la capitale (ou la plus grande ville) de l’état, afin de tester notre fonction,

> carte(US$Frost, + titre="Nombre de jours de gel par an")

Pour le taux d’homicide (qui est la variable par défaut) on a

> carte()

Ça sera notre variable d’intérêt lors de la modélisation de mardi. Enfin, on peut lancer une régression, et représenter spatialement les résidus,

> reg=lm(Murder~.-NOM-USPS-nom,data=US) 
> regs=step(reg) 
> carte(residuals(regs), 
+ titre="Résidus de la régression")

Nous voila équipés pour commencer l’économétrie spatiale…

Pour aller maintenant un peu plus loin dans la modélisation, je vais rajouter une variable qualitative, par exemple l’appartenance politique du gouverneur en 1977 (les données datent de cette époque). Les données ont été extraites de wikipedia, suite aux élections de 1974, 1975, 1976 et 1977,

> GV=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/governor.csv", 
+ header=TRUE,sep=";") 
> etat=strsplit(as.character(GV$State),"-") 
> listeetat=rep(NA,nrow(GV)) 
> for(i in 1:nrow(GV)){ 
+ listeetat[i]=etat[[i]][1] 
+ } 
> indice=which(is.na(listeetat)==FALSE) 
> basegv=data.frame(state=tolower(listeetat[indice]), 
+ party=GV$Party[indice])

On a la visualisation suivante

> library(maps) 
> library(RColorBrewer) 
> couleur=rev(brewer.pal(6, "RdBu")) 
> Z=rep(6,length(basegv$party)) 
> Z[basegv$party=="Democratic"]=1 
> VL0=strsplit(map("state")$names,":")
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,basegv$state)
> niveau=Z[ETAT] 
> map("state", fill = TRUE, col=couleur[niveau])

For those who think more variates should be added to the dataset, some can be found e.g. on http://www.statemaster.com/, like the total executions since 1930, or the date the state joint the U.S.A.

Des calculs en base 12 à la magie des log-rendements

Un des gros avantages dans la culture anglo-saxonne, c’est que culturellement, les anglais ont été de grands amateurs du système duodécimal (en base 12). Par exemple, il y a 12 pouces dans un pied, ou le fait que jusqu’aux années 70, il y avait 12 pences dans un schilling… Et ce que j’apprécie, c’est que ce système pousse à manipuler les fractions (2, 3, 4, 6 étant des diviseurs de 12) alors que le système décimal pousse à utiliser les nombres à virgule. Et compte tenu de la proximité avec les États-Unis, même si ici le système métrique a été adopté, il reste des fonds de ces calculs en base 12.
Par exemple, à l’école à Montréal, les tables d’addition et de multiplication vont jusqu’à 12. Aussi, l’autre soir, je faisais réviser les tables (de multiplication) de 11 et 12 à mon fils, et on a passé du temps sur http://freakonometrics.blog.free.fr/public/perso5/12-18.gif et http://freakonometrics.blog.free.fr/public/perso5/12-19.gif. Car sa première réponse à combien font huit fois douze ? a été ben…. cent ?. J’ai bien entendu protesté (c’est un euphémisme) pour la forme, tout en sachant qu’il est naturel de croire que http://freakonometrics.blog.free.fr/public/perso5/12-01.gif, ou que http://freakonometrics.blog.free.fr/public/perso5/12-02.gif. N’est ce pas ce qu’on pense quand on apprend qu’un prix a augmenté de 20% avant de baisser de 20% le lendemain ? Car c’est la même erreur que l’on commet. Et mine de rien, raisonner avec des rendements, c’est pénible…. (tous ceux qui ont discuté taux d’intérêt avec un banquier le savent).
Pour rappel, on définie le rendement comme la variation relative, i.e.

http://freakonometrics.blog.free.fr/public/perso5/12-03.gif

L’intérêt de cette formule est évident

http://freakonometrics.blog.free.fr/public/perso5/12-04.gif

Et si on regarde ce qui se passe sur deux jours

http://freakonometrics.blog.free.fr/public/perso5/12-05.gif

Aussi, si un jour on a un rendement de http://freakonometrics.blog.free.fr/public/perso5/12-06.gif, et le lendemain un rendement de http://freakonometrics.blog.free.fr/public/perso5/12-07.gif

http://freakonometrics.blog.free.fr/public/perso5/12-08.gif

Autrement dit, si on subit une baisse de 10% puis une hausse de 10% (ou le contraire d’ailleurs), on ne revient pas à la position initiale, car http://freakonometrics.blog.free.fr/public/perso5/12-09.gif(ou ce que je voyais avec mon fils quand on notait que http://freakonometrics.blog.free.fr/public/perso5/12-10.gif). En fait, on perd de l’argent… Damned, c’est pénible ça, non ? D’autant plus que si on répète cette alternance (entre une baisse de 10% puis une hausse de 10%) sur un an, au final, on a perdu 85% de valeur… Autrement dit, avec des rendements d’espérance nulle, et symétriques, la volatilité à elle seule fait perdre de l’argent…
Si on s’intéresse maintenant aux log-rendement,

http://freakonometrics.blog.free.fr/public/perso5/12-11.gif

cette fois, on note que

http://freakonometrics.blog.free.fr/public/perso5/12-12.gif

(on retrouve très explicitement la croissance exponentielle des prix). En faisant des développements limités, i.e. si les rendements sont faibles, alors on observe que http://freakonometrics.blog.free.fr/public/perso5/12-13.gif. En effet

http://freakonometrics.blog.free.fr/public/perso5/12-15.gif

Aussi, si http://freakonometrics.blog.free.fr/public/perso5/12-14.gif vaut 0.05, on parlera (abusivement) de hausse de 5%. Si on regarde une évolution sur deux jours,

http://freakonometrics.blog.free.fr/public/perso5/12-16.gif

Autrement dit, si on subit une baisse de 10% (en log-rendement) puis une hausse de 10% (ou le contraire), on revient à la position initiale (et donc on a une variation nulle) car http://freakonometrics.blog.free.fr/public/perso5/12-20.gif.
C’est probablement pour ça que les financiers utilisent des log-rendements: ils n’ont pas (en toute rigueur) la même interprétation que la variation relative, mais au moins on peut faire les calculs de manière très intuitive.

Exchangeability, credit risk and risk measures

Exchangeability is an extremely concept, since (most of the time) analytical expressions can be derived. But it can also be used to observe some unexpected behaviors, that we will discuss later on with a more general setting. For instance, in a old post, I discussed connexions between correlation and risk measures (using simulations to illustrate, but in the context of exchangeable risk, calculations can be performed more accurately). Consider again the standard credit risk problem, where the quantity of interest is the number of defaults in a portfolio. Consider an homogeneous portfolio of exchangeable risk. The quantity of interest is here

http://freakonometrics.hypotheses.org/files/2016/11/credit-01.gif

or perhaps the quantile function of the sum (since the Value-at-Risk is the standard risk measure). We have seen yesterday that – given the latent factor – http://freakonometrics.hypotheses.org/files/2016/11/exch67.gif (either the company defaults, or not), so that

http://freakonometrics.hypotheses.org/files/2016/11/exch66.gif

i.e. we can derive the (unconditional) distribution of the sum

http://freakonometrics.hypotheses.org/files/2016/11/exch60.gif

so that the probability function of the sum is, assuming that http://freakonometrics.hypotheses.org/files/2016/11/exch76.gif

http://freakonometrics.hypotheses.org/files/2016/11/exch68.gif

Thus, the following code can be used to calculate the quantile function

> proba=function(s,a,m,n){
+ b=a/m-a
+ choose(n,s)*integrate(function(t){t^s*(1-t)^(n-s)*
+ dbeta(t,a,b)},lower=0,upper=1,subdivisions=1000,
+ stop.on.error =  FALSE)$value
+ }
> QUANTILE=function(p=.99,a=2,m=.1,n=500){
+ V=rep(NA,n+1)
+ for(i in 0:n){
+ V[i+1]=proba(i,a,m,n)}
+ V=V/sum(V)
+ return(min(which(cumsum(V)>p))) }

Now observe that since variates are exchangeable, it is possible to calculate explicitly correlations of defaults. Here

http://freakonometrics.hypotheses.org/files/2016/11/exch70.gif

i.e.

http://freakonometrics.hypotheses.org/files/2016/11/exch71.gif

Thus, the correlation between two default indicators is then

http://freakonometrics.hypotheses.org/files/2016/11/exch73.gif

http://freakonometrics.hypotheses.org/files/2016/11/exch75.gif

Under the assumption that the latent factor is beta distributed

http://freakonometrics.hypotheses.org/files/2016/11/exch78.gif

we get

http://freakonometrics.hypotheses.org/files/2016/11/exch80.gif

Thus, as a function of the parameter of the beta distribution (we consider beta distributions with the same mean, i.e. the same margin distributions, so we have only one parameter left, with is simply the correlation of default indicators), it is possible to plot the quantile function,

> PICTURE=function(P){
+ A=seq(.01,2,by=.01)
+ VQ=matrix(NA,length(A),5)
+ for(i in 1:length(A)){
+ VQ[i,1]=QUANTILE(a=A[i],p=.9,m=P)
+ VQ[i,2]=QUANTILE(a=A[i],p=.95,m=P)
+ VQ[i,3]=QUANTILE(a=A[i],p=.975,m=P)
+ VQ[i,4]=QUANTILE(a=A[i],p=.99,m=P)
+ VQ[i,5]=QUANTILE(a=A[i],p=.995,m=P)
+ }
+ plot(A,VQ[,5],type="s",col="red",ylim=
+ c(0,max(VQ)),xlab="",ylab="")
+ lines(A,VQ[,4],type="s",col="blue")
+ lines(A,VQ[,3],type="s",col="black")
+ lines(A,VQ[,2],type="s",col="blue",lty=2)
+ lines(A,VQ[,1],type="s",col="red",lty=2)
+ lines(A,rep(500*P,length(A)),col="grey")
+ legend(3,max(VQ),c("quantile 99.5%","quantile 99%",
+ "quantile 97.5%","quantile 95%","quantile 90%","mean"),
+ col=c("red","blue","black",
+"blue","red","grey"),
+ lty=c(1,1,1,2,2,1),border=n)
+}

e.g. with a (marginal) default probability of 15%,

> PICTURE(.15)

On this graph, we observe that the stronger the correlation (the more on the left), the higher the quantile… Note that the same graph can be plotted with on the X-axis the correlation,


Which is quite intuitive, somehow. But if the marginal probability of default decreases, increasing the correlation might decrease the risk (i.e. the quantile function),

> PICTURE(.05)

(with the modified code to visualize the quantile as a function of the underlying default correlation) or even worse,

> PICTURE(.0075)

And it because all the more counterintuitive that the default probability decreases ! So in the case of a portfolio of non-very risky bond issuers (with high ratings), assuming a very strong correlation will lower risk based capital !

de Finetti’s theorem and exchangeability

This week, we will start to work on multivariate models, and non-independence. The first idea to discuss non-independence will be to use the concept ofexchangeability. A sequence of random variable http://freakonometrics.blog.free.fr/public/perso5/exch-06.gif is said to be exchangeable if for all http://freakonometrics.blog.free.fr/public/perso5/exch-05.gif

http://freakonometrics.blog.free.fr/public/perso5/exch-01.giffor any permutation http://freakonometrics.blog.free.fr/public/perso5/exch-02.gif of http://freakonometrics.blog.free.fr/public/perso5/exch-03.gif. A standard example is the case wherehttp://freakonometrics.blog.free.fr/public/perso5/exch-07.gif, with

http://freakonometrics.blog.free.fr/public/perso5/exch-08.gifand

http://freakonometrics.blog.free.fr/public/perso5/exch-09.gifSince http://freakonometrics.blog.free.fr/public/perso5/exch-19.gif, a necessary condition is that

http://freakonometrics.blog.free.fr/public/perso5/exch-11.gifi.e.

http://freakonometrics.blog.free.fr/public/perso5/exch-12.gif
Since this inequality should hold for all http://freakonometrics.blog.free.fr/public/perso5/exch-05.gif it comes that necessarily http://freakonometrics.blog.free.fr/public/perso5/exch-13.gif.
de Finetti (1931): Let http://freakonometrics.blog.free.fr/public/perso5/exch-06.gif be a sequence of random variables with values in http://freakonometrics.blog.free.fr/public/perso5/exch-14.gifhttp://freakonometrics.blog.free.fr/public/perso5/exch-06.gif is exchangeable if and only if there exists a distribution function http://freakonometrics.blog.free.fr/public/perso5/exch-15.gif on http://freakonometrics.blog.free.fr/public/perso5/exch-16.gif such that

http://freakonometrics.blog.free.fr/public/perso5/credit-04.gifwhere http://freakonometrics.blog.free.fr/public/perso5/exch-20.gif. Note that http://freakonometrics.blog.free.fr/public/perso5/exch-15.gif is the distribution function of random variable

http://freakonometrics.blog.free.fr/public/perso5/exch-22.gifA nice proof of that result can be found in Heath & Sudderth (1995) – see alsoSchervish (1995)Chow & Teicher (1997) or Durrett (2010) and also probably in several bayesian books because that result has a strong interpretation in bayesian inference (as far as I understood, see e.g. Jaynes (1982)).
From the exchangeability condition, for any permutation http://freakonometrics.blog.free.fr/public/perso5/defi02.gif of http://freakonometrics.blog.free.fr/public/perso5/defi03.gif,

http://freakonometrics.blog.free.fr/public/perso5/defi01b.gifthat can be inverted in

http://freakonometrics.blog.free.fr/public/perso5/defi05.gifThe idea is then to extend the size of the vector http://freakonometrics.blog.free.fr/public/perso5/defi09.gif, i.e. for all http://freakonometrics.blog.free.fr/public/perso5/defi07.gif, define

http://freakonometrics.blog.free.fr/public/perso5/defi10.gifso that, if we condition on http://freakonometrics.blog.free.fr/public/perso5/defi11.gif,

http://freakonometrics.blog.free.fr/public/perso5/defi08.gifbut since given the sum of components of http://freakonometrics.blog.free.fr/public/perso5/defi11.gif, all possible rearrangements of the ones among the http://freakonometrics.blog.free.fr/public/perso5/GPD11.gif elements are equally likely, we can write

http://freakonometrics.blog.free.fr/public/perso5/defi15.gifThe first idea is to work on the blue term, and to invocate a theorem of approximation of the hypergeometric distribution http://freakonometrics.blog.free.fr/public/perso5/defi17.gif to a binomial distribution http://freakonometrics.blog.free.fr/public/perso5/defi19.gif, when http://freakonometrics.blog.free.fr/public/perso5/defi50.gif becomes large. Then

http://freakonometrics.blog.free.fr/public/perso5/defi16.gifLet http://freakonometrics.blog.free.fr/public/perso5/defi20.gif and let http://freakonometrics.blog.free.fr/public/perso5/defi21.gif denote the cumulative distribution function of http://freakonometrics.blog.free.fr/public/perso5/defi28.gif.

http://freakonometrics.blog.free.fr/public/perso5/defi33.gifThe idea is then to write the sum as an integral, with respect to that distribution,

http://freakonometrics.blog.free.fr/public/perso5/defi30.gifThe theorem is then obtained since http://freakonometrics.blog.free.fr/public/perso5/defi31.gif, i.e.

http://freakonometrics.blog.free.fr/public/perso5/defi32.gifIn the case of non-binary sequences, there is an extension of the previous result,
Hewitt & Savage (1955): Let http://freakonometrics.blog.free.fr/public/perso5/exch-06.gif be a sequence of random variables with values in http://freakonometrics.blog.free.fr/public/perso5/exch-24.gif.  http://freakonometrics.blog.free.fr/public/perso5/exch-06.gif is exchangeable if and only if there exists a measure http://freakonometrics.blog.free.fr/public/perso5/exch-25.gif on http://freakonometrics.blog.free.fr/public/perso5/exch-26.gif such that

http://freakonometrics.blog.free.fr/public/perso5/exc99.gifwhere http://freakonometrics.blog.free.fr/public/perso5/exch-25.gif is the measure associated to the empirical measure

http://freakonometrics.blog.free.fr/public/perso5/exch-29.gifand

http://freakonometrics.blog.free.fr/public/perso5/exc98.gifFor instance, in the Gaussian case mentioned earlier, if

http://freakonometrics.blog.free.fr/public/perso5/exch-23.gifthen

http://freakonometrics.blog.free.fr/public/perso5/exch-30.gifwhere

http://freakonometrics.blog.free.fr/public/perso5/exch-31.gifi.e. conditionally on http://freakonometrics.blog.free.fr/public/perso5/exch-32.gif, the http://freakonometrics.blog.free.fr/public/perso5/exch-06.gif are conditionally independent, with distribution http://freakonometrics.blog.free.fr/public/perso5/exch-33.gif. The proof can be found in Kingman (1978) and is based on martingale arguments.
Note that in the Gaussian case, http://freakonometrics.blog.free.fr/public/perso5/excccc02.gif where http://freakonometrics.blog.free.fr/public/perso5/exccc03.gif are i.i.d. random variables. To go further on exchangeability and related topics, see Aldous (1985)  (see also here).
This construction can be used in credit risk, to model defaults in an homogeneous portfolio, see e.g. Frey (2001),

 

Assuming a Beta distribution for the latent factor, we can derive the probability distribution of the sum

http://freakonometrics.blog.free.fr/public/perso5/credit-01.gifSince

http://freakonometrics.blog.free.fr/public/perso5/exch61.gifif we assume that – given the latent factor – http://freakonometrics.blog.free.fr/public/perso5/exch67.gif (either the company defaults, or not),

http://freakonometrics.blog.free.fr/public/perso5/exch66.gifi.e.

http://freakonometrics.blog.free.fr/public/perso5/exch63.gifThus, we can derive the (unconditional) distribution of the sum

http://freakonometrics.blog.free.fr/public/perso5/exch60.gifi.e.

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

> proba=function(s,a,m,n){
+ b=a/m-a
+ choose(n,s)*integrate(function(t){t^s*(1-t)^(n-s)*
+ dbeta(t,a,b)},lower=0,upper=1)$value
+ }

Based on that function, it is possible to plot the probability distribution over http://freakonometrics.blog.free.fr/public/perso5/credit-5.gif. In the upper corner is plotted the density of the Beta distribution.

> a=2
> m=.2
+ n=10
+ V=rep(NA,n+1)
+ for(i in 0:n){
+ V[i+1]=proba(i,a,m,n)}
> barplot(V,names.arg=0:10)

http://freakonometrics.blog.free.fr/public/perso5/exchangeable-beta.gif

Those two theorems are extremely close,

De Finetti’s theorem: a random sequence http://freakonometrics.blog.free.fr/public/perso5/dfhs1.gif of http://freakonometrics.blog.free.fr/public/perso5/dfhs4.gif random variables is exchangeable if and only if http://freakonometrics.blog.free.fr/public/perso5/dfhs2.gif‘s are conditionnally independent, conditionnally on some random variable http://freakonometrics.blog.free.fr/public/perso5/dfhs3.gif.

Hewitt-Savage’s theorem: a random sequence http://freakonometrics.blog.free.fr/public/perso5/dfhs1.gif is exchangeable if and only if http://freakonometrics.blog.free.fr/public/perso5/dfhs2.gif‘s are conditionnally independent, conditionnally on some sigma-algebra http://freakonometrics.blog.free.fr/public/perso5/dfhs5.gif

Olshen (1974), proposed an interesting discussion about those theorems, see also in the Encyclopedia of Statistical Science,

http://freakonometrics.free.fr/copecran1.png

The subtle difference between those two theorem is also discussed in Freedman (1965)

http://freakonometrics.free.fr/copecran2.png

MAT8886 from tail estimation to risk measure(s) estimation

This week, we conclude the part on extremes with an application of extreme value theory to risk measures. We have seen last week that, if we assume that above a threshold http://freakonometrics.blog.free.fr/public/perso5/qt01.gif, a Generalized Pareto Distribution will fit nicely, then we can use it to derive an estimator of the quantile function (for percentages such that the quantile is larger than the threshold)

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

It the threshold is http://freakonometrics.blog.free.fr/public/perso5/qt02.gif, i.e. we keep the http://freakonometrics.blog.free.fr/public/perso5/qt04.gif largest observations to fit a GPD, then this estimator can be written

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

The code we wrote last week was the following (here based on log-returns of the SP500 index, and we focus on large losses, i.e. large values of the opposite of log returns, plotted below)

> library(tseries)
> X=get.hist.quote("^GSPC")
> T=time(X)
> D=as.POSIXlt(T)
> Y=X$Close
> R=diff(log(Y))
> D=D[-1]
> X=-R
> plot(D,X)
> library(evir)
> GPD=gpd(X,quantile(X,.975))
> xi=GPD$par.ests[1]
> beta=GPD$par.ests[2]
> u=GPD$threshold
> QpGPD=function(p){
+ u+beta/xi*((100/2.5*(1-p))^(-xi)-1)
+ }
> QpGPD(1-1/250)
97.5%
0.04557386
> QpGPD(1-1/2500)
97.5%
0.08925095

This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days)

> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/250, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI   Estimate   Upper CI
0.04172534 0.04557655 0.05086785

or the decennial one

> gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/2500, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI   Estimate   Upper CI
0.07165395 0.08925558 0.13636620

Note that it is also possible to derive an estimator for another population risk measure (the quantile is simply the so-called Value-at-Risk), the expected shortfall (or Tail Value-at-Risk), i.e.

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

The idea is to write that expression

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

so that we recognize the mean excess function (discussed earlier). Thus, assuming again that above http://freakonometrics.blog.free.fr/public/perso5/qt01.gif (and therefore above that high quantile) a GPD will fit, we can write

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

or equivalently

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

If we substitute estimators to unknown quantities on that expression, we get

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

The code is here

> EpGPD=function(p){
+ u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi)
+ }
> EpGPD(1-1/250)
97.5%
0.06426508
> EpGPD(1-1/2500)
97.5%
0.1215077

An alternative is to use Hill’s approach (used to derive Hill’s estimator). Assume here that http://freakonometrics.blog.free.fr/public/perso5/qt20.gif, where http://freakonometrics.blog.free.fr/public/perso5/qt21.gif is a slowly varying function. Then, for all http://freakonometrics.blog.free.fr/public/perso5/qt23.gif,

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

Since http://freakonometrics.blog.free.fr/public/perso5/qt21.gif is a slowly varying function, it seem natural to assume that this ratio is almost 1 (which is true asymptotically). Thus

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

i.e. if we invert that function, we derive an estimator for the quantile function

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

which can also be written

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

(which is close to the relation we derived using a GPD model). Here the code is

> k=trunc(length(X)*.025)
> Xs=rev(sort(as.numeric(X)))
> xiHill=mean(log(Xs[1:k]))-log(Xs[k+1])
> u=Xs[k+1]
> QpHill=function(p){
+ u+u*((100/2.5*(1-p))^(-xiHill)-1)
+ }

with the following Hill plot

For yearly and decennial events, we have here

> QpHill(1-1/250)
[1] 0.04580548
> QpHill(1-1/2500)
[1] 0.1010204

Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles,

> quantile(X,1-1/250)
99.6%
0.04743929
> quantile(X,1-1/2500)
99.96%
0.09054039

Note that it is also possible to use some functions to derive estimators of those quantities,

> riskmeasures(gpd(X,quantile(X,.975)),1-1/250)
p   quantile      sfall
[1,] 0.996 0.04557655 0.06426859
> riskmeasures(gpd(X,quantile(X,.975)),1-1/2500)
p   quantile     sfall
[1,] 0.9996 0.08925558 0.1215137

(in this application, we have assumed that log-returns were independent and identically distributed… which might be a rather strong assumption).

Delta method and quantile estimation

An alternative (to profile likelihood techniques) to derive confidence intervals is to use the delta method. Consider an estimator such that

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

then for any differentiable transformation

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

The proof of that result is based on Taylor’s expansion (see here or there for more details on the theory, or even on this blog – here, in French or there in English – for some codes in R). This can be used to derive an asymptotic confidence interval for a quantile. Consider the following dataset

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM

It is possible to fit a Generalized Pareto distribution on observations above a given threshold,

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

In that case, if http://freakonometrics.blog.free.fr/public/perso5/GPD10.gif exceed the threshold out of a sample of size http://freakonometrics.blog.free.fr/public/perso5/GPD11.gif, the estimator of the quantile

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

i.e. http://freakonometrics.blog.free.fr/public/perso5/GPD05.gif. Then

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

whilehttp://freakonometrics.blog.free.fr/public/perso5/GPD07.gif

i.e. it is now possible to implement the delta-method to derive the asymptotic variance of the quantile, and also (asymptotic) confidence intervals.

> u=5
> GPD=gpd(X,u)
> theta=GPD$par.ests
> sigma=GPD$varcov
> k=GPD$n.exceed
> n=length(X)
> p=.975
> Q=u+theta[2]/theta[1]*((n*(1-p)/k)^(-theta[1])-1)
> nabla=c(-theta[2]/theta[1]^2*((1-p)^(-theta[1])-1)-
+ theta[2]/theta[1]*(1-p)^(-theta[1]*log(1-p)),
+ 1/theta[1]*((1-p)^(-theta[1])-1))
> variance=t(nabla)%*%sigma%*%nabla

Based on the assumption of normality, it is possible to derive confidence intervals, and to compare them with the one obtained in R,

> c(Q-1.96/sqrt(k)*sqrt(variance),
+   Q+1.96/sqrt(k)*sqrt(variance))
[1] 13.11562 16.82852
+ tailplot(gpd(X,5))
+ gpd.q(tailplot(gpd(X,5)), .975, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI Estimate Upper CI
13.33329 14.97207 17.18094

a short word on profile likelihood

Profile likelihood is an interesting theory to visualize and compute confidence interval for estimators (see e.g. Venzon & Moolgavkar (1988)). As we will use is, we will plot

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

But more generally, it is possible to consider

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

where http://freakonometrics.blog.free.fr/public/perso5/profilik03.gif. Then (under standard suitable conditions)

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

which can be used to derive confidence intervals.

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM
> u=5

The function to draw the profile likelihood for the tail index parameter is then

> Y=X[X>u]-u
> loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
> XIV=(1:300)/100;L=rep(NA,300)
> for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
> plot(XIV,L,type="l")

It is possible to use it that profile likelihood function to derive a confidenceinterval,

> PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
> (OPT=optimize(f=PL,interval=c(0,3)))
$minimum
[1] 0.6315989

$objective
[1] 754.1115
> up=OPT$objective
> abline(h=-up)
> abline(h=-up-qchisq(p=.95,df=1)/2,col="red")
> I=which(L>=-up-qchisq(p=.95,df=1)/2)
> lines(XIV[I],rep(-up-qchisq(p=.95,df=1)/2,length(I)),
+ lwd=5,col="red")
> abline(v=range(XIV[I]),lty=2,col="red")

This is done with the following code

> library(ismev)
> gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)